{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--BOOK_INFORMATION-->\n",
    "<img align=\"left\" style=\"padding-right:10px;\" src=\"figures/PHydro-cover-small.png\">\n",
    "*This is the Jupyter notebook version of the [Python in Hydrology](http://www.greenteapress.com/pythonhydro/pythonhydro.html) by Sat Kumar Tomer.*\n",
    "*Source code is available at [code.google.com](https://code.google.com/archive/p/python-in-hydrology/source).*\n",
    "\n",
    "*The book is available under the [GNU Free Documentation License](http://www.gnu.org/copyleft/fdl.html). If you have comments, corrections or suggestions, please send email to satkumartomer@gmail.com.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!--NAVIGATION-->\n",
    "< [Interpolation](05.09-Interpolation.ipynb) | [Contents](Index.ipynb) | [Uncertainty Intervals](05.11-Uncertainty-Intervals.ipynb) >"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.10 自相关\n",
    "\n",
    "自相关是信号与自身的相关性。这告诉信号如何在时间或空间上相关。这可以用来查看信号的周期性。为了证明这点，我们将首先使用周期为$4\\pi$和幅度为2的正弦信号。图5.17显示了信号在上面板，自相关在下面板。使用`plt.acorr`函数绘制自相关。我们使用`plt.grid`函数在图中显示了网格。使用`plt.axhline`绘制了网格水平线0和$e^{-1}$。自相关显示了一个具有周期为$4\\pi$的很好的周期性行为。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEKCAYAAADJvIhZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xl8XXWd8PHPN/veJE2T7k032gKyFiiKEoEqm6LPoKKjg8z4MM44wqgzDiiugwuPjI7bOCIuqDwgD4IwgkApDbJIpRsUutDaNWmaNM2eNPv3+eOcm9yk9yb35p577s293/fr1VfvPffk3N9Z7vme3y6qijHGGOOHjEQnwBhjTPqwoGOMMcY3FnSMMcb4xoKOMcYY31jQMcYY4xsLOsYYY3xjQccYY4xvLOgYY4zxjQUdY4wxvslKdAISpaKiQqurq6f0t93d3RQWFnqboCSXjvsM6bnf6bjPkJ77He0+b968uVlVZ8XynWkbdKqrq9m0adOU/ra2tpaamhpvE5Tk0nGfIT33Ox33GdJzv6PdZxE5GOt3pm3QiYfVt6+juav/pOUVRTlsum1tAlJkjDHJxep0PBQq4Ey03Bhj0o0FHWOMMb6xoGOMMcY3FnSMMcb4xhoSGGNiZo1oTKQsp+OhiqKcqJYbkyqsEY2JlOV0PLTptrVsOdTK//qvFwF47rNvZ0F5QYJTZYwxycNyOh5r6ugded3SbU95xhgTzHI6HjvaHhR0eizoGFN9y2OA1e8YhwUdjzV29o28bnVzOlbJaozV7xiHFa95rLGjl+JcJ5YHitesktWkOmssYyJlOR2PNXb0srSyiO317bRa8ZrxUSJz1Os/U8OZX3mKf3nHKdz51Btx/S4zvVlOx2ONHX3MmZFHWUEOLd0DiU6OSSOJzFG/WtcGwFkLyuL+XWZ6s5yOxxo7erloWQXlhdkjdTomMlb3NX1tO+QEnTMWzEhwSkKzayt5WE7HQz39g3T2DlJZkuvkdKx4LSpW9zV9bTvcxtJZhZTkZSdlJ2m7tpKH5XQ81NjhtFybXZJHeWEOe5u6AOfHFu4py0/J+LQXLk1m+lBVth1uo2ZFJcDItfSH7Q38w71b+P0nL+L0ecmZAzL+s6DjoUa3Y2hVSR5lhTkjDQk23baW9hMDnPmVpwD47nVncc1Z83xPXzI+7VnAmf7qWk9wvLufsxaWjlm+vKoIgD1NnRZ0zAgLOh4aDTq5lBfk0NozwPCwkpEhI58BHAvqy2Omj2TMKY5PRyJy1NsOO/U5Zy8YG3QWzSwkK0NGcvwmvGS/tryU1HU6IrJARDaIyE4ReV1Ebg6xTo2ItIvINvffFxORVhib0yktyGZoWOnsHQTGjlTQZEFnWkrGnGKwu68/D4ClswoB2PXvl3Pgm1fF/aa17XAbuVkZrJhdPGZ5dmYG1RWF7Gm0oDOZZL+2vJTsOZ1B4DOqukVEioHNIrJOVXeMW+85Vb06Aekbo7Gjj4KcTIpysygvdJ4uW3r6mVGQPRKQsjIk7XM60dbjWMfDyKzf2UiGwLXnLuCOJ3ZxvLufeaX5cfu+8edx+ef/AIx9Ol82q4g3GjvjloZIJUu96lSk2jBCSR10VLUBaHBfd4rITmAeMD7oJIWjHb1UleQhIpQFgk53P4srCkeCzilVxTR19k60mZQXTcC5+29Wc9mpVXFMTepYv7OJ1YvKWVbp1KU0d/bFNehE8nS+vKqIdTsb6RscIjcrM25pmcym29by4bs38vzeZgpzMnntK+9ERBKWnqlIlVxPUgedYCJSDZwNbAzx8YUi8gpwBPgXVX3dx6SNaOropaokF4DyAifoBPrqHO3opbQgmwXl+ew71u1bmiLJVSTj015FUQ7DCr/bVm9BZwKhzu+fD7QA0NyV+Bz1ssoihoaVA809JxW/+e1gSzci0N0/xLGuPiqL8xKannQ1LYKOiBQBvwX+WVU7xn28BVikql0iciXwO2B5mO3cCNwIUFVVRW1t7ZTS09XVFfJvDzT2sLQ0g9raWo71DAPwpy2vktWUzev7einKGGag8zhHWgen/N3RmijgXL0km2tPcQLOZOkJt89e+8XlTn3ETc9009EPv3+1gd+/+tjI5yU58L1LCuOejoBI9zsexyZwDMYLPgYTnd/nN79KZmN21N/rxbkO/H1bxxAAj9Ru5PzZibvdDA4r9a0nqC7OYH/HMA+ve4EV5WNzXn5d47HwOn2J2OekDzoiko0TcO5V1YfGfx4chFT1cRH5LxGpUNXmEOveBdwFsHr1aq2pqZlSmmpraxn/t6pK+9NPcObyRdTUrKKrb5B//eOTVC5YQs3FS/mP7c+zdG4OZy8q45lDb3DhRW/1p7jhicfCftSbW0ZNzXkRbSbUPscjTYHv6AizTkc/3qUjAsH7XfF8+BZG8UhTRMdggmNZMW8xNTXLov7eiM91BOexd2CIL//pCXIqFlJTc0rUafHKgeZuhp+q5erVS/j+M3spW3AKNectGLOOp9d4lMqfeyqiYbO8Tl8i9jnZW68J8FNgp6p+O8w6s931EJHzcfbpuH+pdLSfGKB/cJjKEifLXpiTSU5mxsioBEc7epldksesYqf4LRkaE+w6mvgK3ulk021reevyCorznGe173zgTF9ah01VMhSv5WVnsqCsIOHNpg+19ABw4ZKZZGUI+4/7V8Qdic++cyUAT3/6bQlOSfwlddAB3gJ8BLgkqEn0lSLycRH5uLvOtcBrbp3O94DrVFX9TmhgNIJAnY7TmMAZf21gaJjmrj6qZuRR6X6eDM2m61pP0Nnr/6CkyThMSiRUlR1HOrh0ZSUZAvt9rJubinhXPEd6HpdVFiU86Bx0g87iWYUsLC/gYJIFnef2NDNnRh5LZxVN299HpJK6eE1VnwcmbGKiqj8AfuBPisI76rZOm10yWjkZGGm6uasPVeezQOVlU0figw7AG42dnLuo3Nfv3Pi5yzjzK0/xnrPncvt73uTrd8fiWGcfx7v7OXNBKVsOtbH/eE+ikzSh5jg/2Gy6bS0337+VLYdaee6zl4Rdb3llEc/vbWZwaJiszMQ85x463k1OVgZVxXlUVxSyvzl5zt3QsPL83mbeeVoVIjKSc65r7eGiOzbw79ecxkcurE5sIj2U1EFnOgnuGBpQXphDW0//SMfQ2TNyqQwUr/lU9BGuf0J5YTYt3QPsOup/0NnT1ElX3yDnLEyeYfDDtfIryYFXa5zXrzc41YenzZ3h3rgS9/Qe6LsRTk5mhi/Fa81dfVQU5Yb9PPi4LnP78YD/fU4OHu9hYXkBGRlC9cxCXtp3HFVNimbT2+vbaT8xwEXLZ41ZPq80n1nFuWw91MZHLkxQ4uLAgo5HmtygE6izASgrzGFnQ8eYgFRemIMIHOvwp6/OptvW8pGfbqSlu5/HbnrryHJV5YwvP8XuBNTrbDnoDJsyUdDxuzNfuKKo4NZjO444QWflnGKWVBSy5WBrXG9c4Y7BRA588yoAPv/wdh7f3hCPZI3R3NnPwpkF4T9Pkp72h1p6WFTupLO6ooCe/iGaOvvGPCT6KdRDzk33beWr//P6SDAWEc5eUMpWd5ihVGFBJ0bjL56VX3gCcG4YV5w+h9buoJxOSR5ZmRnMLMz1rU5HVXn9SAeXraocs1xEOGV2Mbsa/A86Ww+1Ul6Yw6IJblaBH97B491c/K1a/s+1Z/D+1QvCrh9PgVxFTlYGC8sLKMnLpnpmAV19g3Ht77HptrVc/K0NNHf20d0/xCOfeAtnLiidNJcDUFGUS2vPQNyLtJq7+jhnUfLkWENRVQ619HDh0pkAVM90mpvvb+5OWNCJNBifvbCMp3Y00trdP9LhfLqzoDMFNz3THbY5a0Bzl3ORtJ0YoKG9l+xMocztMFpZ7F/QaWjvpaW7P+QovytnF/M/rxzxvZhhy6FWzl5QGtF3zpmRj4jT6CHR+geHOXVOCQCLZzm9/g8093Dld5+Ly2CN/YPDHG7p4eoz5vLoK0fY2dDBmeMG1Qynws1xt3T3j7So9Nrg0DAtPf3MSvIK7uaufnr6h1jo5nQWVzhB50BzN2uWzExk0iZ1tjty97bDbbx9ZeUka08PFnSmIFSHvVDKC7JRdZomVxbnkZHh3GQrS3J9azL9+pFAPUTJSZ+tnF3MvRsHOdrRy5wZ8RsuJVh7zwB/OdbN/zpnfkTr52RlMKckj7rW5Kj4PdU9jotHnpa74laEdLi1h2GFt50yi/U7G6Nq4l7hPhUf6+qLW9Bp6elHdTTARcuvMcUOtTgt1QI567ml+WRnCgfi1BDEyxGjz5g/gwxxSgcs6JhJBbLDOxs6WFA+WpQ0qyiXnQ3jB1aIj9fq2xGBVXPGBp3gH8aF33hmZHm8bwBbD7cCJw+DP5H5ZQXUtSQ+pwOM5HTmlTk3rn3N8Wt6G2iSvWRWIStmF7MjimsmEAjiWXfS3Olse6KGBBFtJ871Owfd4LKw3HlQyMwQFpQXcCBO587Lh5CCnCxWzi5JqXodCzpxFBhpuqmzj/OqR1uIVZbk0tzVz9CwkpkR32Kt1490sKSikIKcsafa7wre8U9/H7rbGUIvkiA3vyyfl/bFt79vpJX2H/vlppH1F8bxxgVwwO1LsqSikFVzSnjULQqNpJFFIBDEs9l0oHXcREFnKo0hvHaopQcRWFA+mptfPLNw5Pgms+DfTXBd3nQecdqCThwF6nBgbFPqyuI8hoaVlu7+Ma3d4uH1I+2cvzi6JtHxKPaIJcjNLy/g6LZ6+geHycmKT6X4ptvW8tkHX+HZN46NdPSdSHNXP2ctKGN/HIPOvuZuSguyKS3IYdWcEu7deIj6thN897qz+eu7N/Kzj67mkpWhB0MNBKB4NpseDTrh63SCr59IGkDEw6HjPcwpyRsz7FR1RSEv/KV5ZJJFvwSOgQCherCPP5bJ0vrPSxZ04qSiKGckpwNOH52AQKBp6uyNKehMVnZ8vKuPhvbekPU5kUiWC3t+WT7DCg3tJ1g0M36DfR5p62VuaT5DwxrRvi+uKOCPe47FLT0HmrtHKr1XzXFGaN7Z0MnWQ61kZsiY3PN4RblZ5GbFt6/OSNDx4MEpVEDy6qHnYEvPmOLt4N/Nks89PrI8uE9WvCnOOXrHqVV8+wNn+fOlScKCzhSU5IRuTDD+R3Kif2jk9dicjjfjr032FBRoRHD63Ok9P/2CMueGUdca36BT33aCU+eW8PA/vmVk2URP54sriugfHGZGfjbtJ04eTijWPkX7m7u50G1dtWK28+Cws6GDF/5ynDPnz6A4L/wI0iJCRVEux+NZp9PVT05WBsW58bmNxPrQM/6hbLKcVqQNhLzS1TfIQ1vr+eOeY9O2qGwqLOhMwfcuKeTXBws5eLyHdZ++OOx6+TmZ5GdncmJgaMzwOCND4cSxvD34BxZN/Ukyml/mlMUfbolfCzZVpb7tBGujmLunusIJhqdUFbHpYCv/8b4z+fQDr3DP357PxafMmuSvJ3aif4iG9t6RnE5RbhaLZhbw5/0tbK9r4xNvn3z06Iri3JNGvvCyZVVzVx+zinIjbm7vd/1OonLq0e5nspQo+MWCzhS9Vt8x0tkslPE/7g/c9RLgXJCBcar8Hmk6OD3JUMEbqTkz8sjMkLj21Tne3U//4DBzZ0TevPhDP3GC+csHnBZ5n37gFQAa22MfbSJQyV1dMZqzWzW7hCdePwrAm5dWTLqNisIcjoxLi5d1BM1d/VHl5oKD2iPb6rn5/m08+PELufa//xT1dyezwH7e8PM/s2F3/IpfpysLOlPQ0acc7Zi4rmSiH3d+TibFuVkjQ+ckgt8VvLEMa5OVmcGcGfHtq3OkzQloc8dN7xyuKHUiRz04r4FWcYGczviHmA/+ZPQhJlwOpaIol1fr22NOSzjNnX3MjiJIB7t0VRU5WRk85sNQPYnS1TcY8zb8Hg7KDxZ0puCgOxviqVOsoAeYVXJy0Uei+HFhb7ptLT94Zg93PvUGu/79cvKyo5vAbn5ZPofjmNMJF3S+d0lhyEmuJgrUXgSdQP+fQE5nKjmUiuIcWrr749ZCq7mrj9PnTe03UPOtDfQPDvPzFw54m6gk0tk7SE5mBv1Dw1PeRuCB4j+e2s0PN+zljduvSNhI3V6xoDMFBzudi+i0KVbQB25Y+451x9Q82asisuDv/fRvtvHHPcd4+fOXeT40TlNnHyV5WVEHHHAaE8SzpVh9mxMo5pXGPjKDFznYA83dVBbnUhRDJX1FUS5Dw0prTz8zY+zAOd7wsHK8u3/KHUMjuW79fpovCfo6L+q+uvsHueqMOXzHbZ0WS4nC3FKnBWdjZ58n12gixSXoiEiWqsaet0xSBzuGWVCez4z86OefD2eiH+FEP4DSgmyuOH0O9/35kCfpOH9xOQ9trWdfczdL3fHFvNLUMfUhWeaXFdDY0Uff4FBcpvk+0naC/OxMSgtiP6de5HT2N3ePqc+ZikBAON7tfdBpOzHA0LDGPBpBKP/6zhURNZSYzEQ5+EDguOW3r/L7VxvY/IXL+NPzz42sM1HOMtIHxa7ewTEPDbGUKARy4EfaTqRX0BGRu4CbVDXsr0pEFgP3AWtiTFtge5cD3wUygbtV9ZvjPs8FfgmcizNN9QdU9YAX3x3OoY5hzlniXzPkyYpWAn04JhPJxR3oSPrn/S3eB53O3pGZVaMV6E1e33qCJR6nC5wf89zSPE9yd0fbp15sGm0z34kEj0pwSlWxu8ybolQv++iMt6cx8jHmJnog+/PnLuPULz3BX1+wiC9cfWrIv3/n6bO5/+XDvLj3+MSzRYYw0YOiqtLVN0hR3ugtNpaWo/NKnYe1QDHwdBZtTudjwBoReb+q7hr/oYhcC/wEmHplx9jtZQI/BNYCdcDLIvKoqu4IWu3vgFZVXSYi1wF3AB/w4vtD6egdoLFHJ+1w6WfrsPHjqo0XmGMlEosrCqkoyuXP+1v44PkLY03aGOOHA4rG/KC+OvELOpE/QYY7vwU5mRzv7mNgaJjsKZS9e3nNBAJJcN3h7e85nY//egs/v+E87nxyN+WFOfzq7y6IPp2dk49GMFV7opjaeqIHsvq2E/QODLOsMvz18q//z2lxeMMvXnYWTDJ6fKT6BocZGNKYikeDBQbkrQ8RdLxsBu+HaI/I14BbgU0i8klV/TmAiOQA/wn8PdAKvNej9J0P7FXVfe733A9cAwQHnWuAL7uvHwR+ICKiqqFGmYjZzsCozSGmCggW7mTHo6XYitmR5XQiISJcsLicP+9v8Wyb4Dz5NXX0jXSMjVagr84/3rslZKugWH9g9W29UTUMCfdd9/35ELc+tJ0mn8reJ7rpB3I6n3toOzffv23MZzf8/GXysjI4kpM5paktAoFsVhyK1/Y2dXkyLuFeN3hNFHTi9WDY7V6jXgWdwtwsSguyQ+Z0pttQOVEdEVX9gojUAr8G7haRS3CKvu4GzgBeAD6oqnUepW8ecDjofR0w/rFsZB1VHRSRdmAm0OxRGoCTnyZu+LnzZOTl08Rk4zKFMq80n5K8bE9boJ1XXcZj2xuoa+0ZyWHEqv3EAP1Dw1Ma9if42IdrhhrLD6x3YIjmrj7mejC9Q6ATcGNHr+dBJ5oca/Ax6w4aGSNY7+AwvYPD1LediPo8B7Y91TqdcNdrUW4WXX2D1LeemHBG0kiMBJ045Iwn0+Vx0AGYOyOfI22J62bhlaiPiKquF5EzgV8BH3L/DQO3A19W1am3DzxZqEed8ffjSNZxVhS5EbgRoKqqitra2ogTMtHTRDTbiaTfRzRZtFnZzvffeVE2ELoSPJr03fRM90j6Lrpjw8jy4mzl+0S+nfHq3RZ/zXX7qK2NrtFDpAElmv0M1tjtpK3j6AFqa+vHfNbV1RXVdg+5zenXv7iZjn3ettPx4noN5TdPvci5VaNpjWSft+zuJ1Ng659fIGMK9WDhrte9rUPcvnGQh9a/yFmVsR2/517dQ3EOvPLyizFtZyLhjlOgW8WBvbuo7dzryXflDvXyRl1nVNfBZOtGe317YapntQs4xugNvx34o8cBB5ycTfAcxfOBI2HWqRORLGAGELJsSFXvAu4CWL16tYbqfxHWBGW90WwneEBBL4raLjp9MTU1K2LeTkC4GVE7B4SPPuH0HZlK7u75Pc3wwkbefsHZXBDtbI0RlrNHdT6DvLC3GZ7byCVrzjlplIna2tqotnu8q48vvfg0FQuWUvOWxdEnxqPrLNJjlpkhaOmCMddQJPv82LFXqDh+jEve/vbI0xSBs08McPvGp8irWkzNxUsn/4MJ9rM7s5hV84Samgun9PeTqSjKCXucXtp3HF58iQvPPYs3L5t89IhIbGh/jYe21p/8nTFcM9Fe316IOui4uZzfAMuBJ4GHgTuBJ0TkDuALHgafl4Hlbou4euA6nJxVsEeB64E/AdcCz8SrPicZzCvNZ9WcYj69dgVXfu85VkbYcs1LUynKaup0igXiNYtlLAKVs14Uh5UX5pCTmRFxs+lwlcDjxavPyvLKIrZPYdSC5q6+uDSXnpGfTVVJLnsaI2tMMFGx8t6mLq46Y44n6Trwzavo6hvk9C89yb+84xT+6ZLlE64fqNMp9LJ4rTSfzt5BOnoHKJlgsNdkF1XzGhH5BM7NfQnwOVW9ws09nAu8CtwCPCcinjR7cvv6/BNOcNsJPKCqr4vIV0Xk3e5qPwVmishe4NNuGlLW2lOr+OOeZjYfcsb7Wjnbk4aCcRcY3HSqDQni6UjbCUSgakbsaRMRKktyIx5/baKAc+CbV438i1crpNPnzeC1+naifU5zxl2Lz7lcXlnM3qbImk1vum0tP/7IuQDcesVKAH7yN6v5w81vo/3EwKT1OdEE86LcLOaX5Uc0bfhInU6et0EHoGFcvU64fUjWoXKiPSLfBw7hNBYYGaVPVfeIyBrgP4BPAFtxKvNjpqqPA4+PW/bFoNe9wPu8+K7p4J2nzeYXLx7g7uf2kZuVQXWMla1+aezopTAn09Mnv2Cx/MCOtJ1gVlGuZ51Oq0ryIpoILtEqinJ407wZPLi5jqMdvSPNciPR3DXa98dryyqLeGDT4TGt6iZqFvyh8xeSmSFcd/5Cvr3uDV7ad3ykAn+ilmswthXiGV98LOyUJQErZxezO4Kg09nrBB0vp30I7iAa3GJ1021r+asfvciuhg66+4f46fWruXRV5KOl+y3aI/II8Leq2jr+A1XtBz4pIutxch8pJR7jk4Xb5kSzCn7yvi3A6Lzvyz7/h5HPkqFNfribQ05WBvOnWHwV7jjlZ2cyMDTM9i+/k/ycqQeMwORtXpldksfOhg7PtjcVkfTGB9h80Pkpb69rjzjoqCrHu/qpKI5TkV9VET39QxxpH20BOFFDnu317SybVcSM/GzOWVjGS/uOj4zmMFnQCRZunL1gK2YXs2H3sUlHxohHTidwLMb31ekfHGZ7fTvvOmMuv91SF9fZbL0QbZPpSfvfqOrvRGTz1JOUnIJ/qF5VvgVv88ZfbmLvsS6e+Yyz3caOXi74+nq+8u7TuP7N1SPrhWt84FWb/Fg7tYb72/7BqTWXhvB9YjbsauKGX7zMtsNtE04zEU643v+xBvCqkjw27G6aUv8Xr0SS/uD9v/FX436yT5x8LMYfrx8/u48fP7vP8wee5ZXOU/yexs6I6tleO9LB25Y78xetWTKT/1z/BlsPtlKYk8mcKY6CHc6K2SUMDSv7jnVP2Cm7u2+QDHEejLwyqziXrAw5qa/OzoYO+geHuWRlJet2HB15IE1WcSnrUNXDk69lgq2aU8K6nY309A9SkJPFDvdJeaWHHT8jEXzzmKhIYyq8bkRwzsIyRODlAy1TCjrx6lQ3e0YuPf1DdPYNJnWFbyT7GbyOH50Qg6+5j7p94SZzrHN0tOsLlpSjT8Nj2xtYMbvY86C/wi1S3H20c8Kg09k7SGFulqffn5khzJ6Rd1LQ2erW756zqJTqisKRuZiSlY0ynSRWzSlB1bmYz15Yxq4Gp9x45SRD3MRTcAB6151PMJxTyGM3vfWk9SJtgeV1I4IZBdmsqCrm5QPejp4Qq8DU5I3tvZMGnVScLyUWUw1gp8+bMeY67Bsc5tW6dqpveczTnNiSWYVkZ8qkjQm6+gbjMo333NKTO4huPdzG7JI85szIp3pmIVsPn1T7kVQs6CSJU93gsrPBCTo7GzqYV+rtSNaxmFsoPHukK+TcLJHeKOLRcu286nIe2lLH4NBw0swzMhJ0OvpYPkll+6bb1vLP92/luT3NvPz5y+Iy702qE3F+P37kxLIzM1g6q4jdRyeus+vqHfS0PidgXmn+SUNUbT3UxtkLSwGonlnA7189Qv/gMG/+5vqkHJPNgk6SmF+WT1Fu1kgF9K6jHRGPHu2HuUUZ9A4MUt92ggXlU2sxVxWHPjqrq8v41UsH2XW0k9NDjIfndRFhJAJD4UTSV2d4WHl+bzMXLa+wgDMFOZkZLCjPj1uryFBWzC7m5UnGJuzuH4xLmuaW5nG0o3dkbLrmrj4OtfTw4TVOL5XqikKGFQ639iTtmGwWdJJERoawak4xOxs66B0Y4i/HunnHqbNPWi9RxTFzi5xcxJ6mzikHHa9zOsEB5ervPz+yPPhJLhE/vMAUzo0RBJ0dDR00d/WPVISb8HZ89Z2c+sUnKcjJpMcdT65/aJi/BE2G6IcVs4t5ZNuRCTtpdvYOUhyHnM7c0nyGhpXzvvY0Ld2j1/DXH9/F1x/fRalbMnIgiVuwWdBJIqvmlPDQlnr2NDqj7IYabSBR2eI5hU7Q2dvUxSUrI+8DcNMly1haWcTN92+jcopz6YQTr4ASSwAPDoTfenI333py98g2Q527wGyob13uzVAp0YikpWLwsUh0/VOgAr0nzACmfgk0JnjjaCerw0zV0dU3yNxSb3P2wddWcMAJ1nZiAIADSdyCzYJOElk1p4SuvoOs29k48j5ZFOUIFUWRD08S8NK+Fkrcp69Zxck3BE6waEZxDieSQBiqyO/8r6/3vaw91HepKqtue5wPrqnmS+867aT1//HezWw+2MpLt14al+bgEwW2+iQYYTn43F373yP94086d+NnDfVCpA9TxXlZk+Z0AjnDkpyx40H6wYJOEgkEmYe31pGXnUH1zNimK/Zlr1T1AAAfAElEQVTa8sqikBNsTdR5c9vhNlbMLiY3K4OSOBQ3TEfJWtYOzjA+M/OF+taT523pHxzmuTeaueqMOXHrfxS4cbd293P2v6/jtqtW8bG3LgGY8pTsXubEIj133X3xqdOJxOIomk1PNuJ9PNhdIImsqComQ+BwywnOnD8j5kmsvLassojfba0/qdPjptvWsvr2dVy6soo7rj1jZHmg8+ZTO45SVeLNVNAm/mbmZ4ScoXLTgRY6+wa5ZGVl3NNQWpBNcW4WdUHBr771BJkZwtBw+HHivMitxmp4WOnqj0+T6UgsmlnItsOtvs5eHI3kaGNqWH37OlZ98QkCv6dX3D4Gq29fl9iEBVleVURn3+BJ44q1nxiguaufxbPG5sxWV5eRIU7T4UQN9Jmu/V1iUZEnI0Fn9e3rqL7lMapveYwP3b0RcEYviPd1KSLMLy/gUMto3cSRthMjLQOTWc/AEKreDoETjcUzC6hvPcG9H1sDwJffdWpC0hGOBZ0kkcxFLgGBEXv3jitiC4z1tKRibNB5+521I0F008HWkZuXVzesSEbX3XTbWvZ/40oKczL56Jurp92IvIkwM19o6xmgq28wodflwvL8MUGnvu0Ec0vzkv4cjk5V7X8fu4qiHBbNdJpN//yF/QBcdmpyDf5pxWsmYsuqnKCzp6mTi4JaW+075gShJeOGkY/3DSu44nbroVbe+18v8t3rzuKas+aNWe94dz/d/UMsmlkQ94r6RLfw8sLMfOdZNFS9jp8WlhdQu/vYSHFufdsJzl1Uxv/7+JsTmq7JBEaYLsz1btw1iH4Q14e21rNqTolnU857xYKOidisolxK8rJOyunsO9ZNZoawcIr9d7xw5vxSZhbm8MyuppOCTmAARD/SF/zj7+4bZM3X13Ppqkr+87qzR5Yne2CqyHPq3urbEtvsdmF5AX2Dwxzr7GNmUS5Hg0adTpRIzl1ghGmv++lEO4hr/+AwOxs6qL7lsbAj15ck4JKzoGMictMz3XQ84UxrdO/GQ9y70WlJVFGUwwWLZ7KgLJ+crMSV1mZkCDUrKnl6Z+NJQ+IcdotoFvk899DF39pAZ98gv9t2hN9tG51lvaIohx9/5Fz+/leb+e0/vJlzF5X5mq7JzMx3g06Cczrz3YeEQy09DCsMDqunU1BMRfCNf8eRDq783nN8/b1v4kMXjM5b2dWbuOK1cKUISuhGFrW1tfFNUAhJG3RE5FvAu4B+4C/ADaraFmK9A0AnMAQMqupqP9OZLsI1rWzu6ucvx7pOKlrzW/ATXmCOIXBu8B9ZUw3gezHDRMWLr9W3kyGjY+4lkxm5Qk5mBnUhWrD5aWFQ0Am0fEx0TifYqjnFVM8s4PHtDWODzkidTtLeXhMqmY/KOuBWVR0UkTuAW4F/C7Pu21W12b+keS/Zi1wmcuB4Nxct879HfbCJbvAHW7qZXZJHnodzm8Rqe307yyuLY5p8Ll4yRJhTmkd964mEXpfzSvMRtwtBIOea6JxOsPO+9jTNXf0cON4zZhieQLCxoBNa0h4VVX0q6O1LwLWJSosfkmHWz6nqHRg+qbk0JE8gPdzSw8Ikm9b7tfp2Lj4l/v1dpmpeaT71bSfYdNtaNh9s4a9+9Cfu+si5vOO0k8cDjJe87Exml+RxqKWH3OxA0EmeJtPhHnTiMWtoKpkuR+Vvgd+E+UyBp0REgR+r6l3+JcsELKk4uXgtWQLpweM9XHxKcg2o2dzVz5vmJV/RWsC80nyefcMZF+7VunYAzphf6ns6FpQVcLilh8LcTErysihO4knxxvO69VqqSGjQEZGngVCPTp9X1UfcdT4PDAL3htnMW1T1iIhUAutEZJeq/jHM990I3AhQVVU15Uq0rq6uhFTAJbPGPa9Qezg5u301dfYx1NE4pXMWz3M90LSP2tqDcdl2LLq6uhhs76epc4B1z2zgqdf6KM0Vdm19iV0+pyW7v48dx4fo725nRrbG9Xfn5bnOyoA/Pf+cJ9uKRklO6PrXkpzQjQYScS9LaNBR1csm+lxErgeuBi5V1ZBjX6jqEff/JhF5GDgfCBl03FzQXQCrV6/WmpqaKaW7traWqf7tdFXyzGMhL+a87AyyMjK45p1vT+wwN09MPLT9xeeeRs24ptSRiOVcVzwfei6fgpxMegeG+PBVNUlZp1NbW8tbzl7Gw3tfYfkZ59O49WVWLymkpuY839PyyuAeXlz/BmXksWJ+fNMQ9bme4Joryc9JyD0i2sE7E3EvS9riNRG5HKfhwMWqGrLDgIgUAhmq2um+fgfwVR+TmTa+d0nhmIvzr370It19g8wqzqWtZyDh46qFqz8qycuio3eQRQkYPDW4eHFoWLnwG+s5Y34pw6rUtfYkZcAJCLQS23W0k33HunnPFAK2FxbOzEfV6QuW6MYq0bBGBOEl85H5AZCLU2QG8JKqflxE5gJ3q+qVQBXwsPt5FvB/VfWJRCU4nbz7zLl86dHXOXi8h3eclvhhNoJv8E+8dpSP/3oz99+4htfq27n9sZ0J7bgKcMHXnZZOT7vTVoAzvHyipw4OZ36ZE3Seev0oAGfMP3lWVj8sCGrmnkwt1yD8g05OpljQmUDSHhlVXRZm+RHgSvf1PuBMP9NlxvaJOTEwxCPbjvDItiNJcwN9y7KZZGUIz75xjJ4+Z7TfsoLEVkBPh7H1gs2ekUeGMDK3UyIaEcDYUSSSLegErnVV5fQvPcn7Vi/gy+8+jevu+hMTDISd9pI26Jjklew30OK8bM5dVMazu49RWZLLgvKChBf/TTfZmRlUleTR0N7L/LJ8ygv97y82frK7m+7byk33bU2ah5sAEWFpZRF/cccg7OobpDLJJyxMpORsbmRMjGpWVLKjoYNX69p9H/5murvpmW6qb3mMhnZnps661hMJmWYj2R9ugi2dVcS+Y85o6/GYNTSV2JExKWf8XPJ/eO1oUtefJJuJhjwyoS2dVcjDW+vp7hukq2/QOoZOwHI6JuVMpydkkxqWumMP7m/upqsvcbOGTgcWdIzxQbJPPGZis7TSCTq7jnbSOzBsxWsTsCNjopYsY6pNJ1asl9oWzSwgQ+DVOmcg/EILOmHZkTFRsxuo8cN0erjJzcpkYXkBrxx2go7V6YRnR8YYM0a48bv8vtlPt4ebpbOKeG6PM8OK1emEZ0fGpJzp9IScjMYPeWQis7SyiPW7mgDL6UzEjoxJOdPtCdmkhqVBc0pZnU541nrNGGM8sDRoynYrXgvPgo4xxnggOOhY8Vp4dmSMMSZG48eJu/AbzwDYKBghWE7HGGNiZKNgRM6CjjHGGN9Y0DHGGOMbCzrGGGN8Y0HHGGOMbyzoGGNMjGwU8ciJanpO5i0ix4CDU/zzCqDZw+RMB+m4z5Ce+52O+wzpud/R7vMiVZ0VyxembdCJhYhsUtXViU6Hn9JxnyE99zsd9xnSc78Tsc9WvGaMMcY3FnSMMcb4xoLO1NyV6AQkQDruM6TnfqfjPkN67rfv+2x1OsYYY3xjOR1jjDG+saBjjDHGNxZ0oiAil4vIbhHZKyK3JDo98SIiC0Rkg4jsFJHXReRmd3m5iKwTkT3u/2WJTqvXRCRTRLaKyO/d94tFZKO7z78RkZTr7ScipSLyoIjscs/5hal+rkXkU+61/ZqI3Ccieal4rkXkZyLSJCKvBS0LeW7F8T33/vaqiJwTjzRZ0ImQiGQCPwSuAE4FPigipyY2VXEzCHxGVVcBa4BPuPt6C7BeVZcD6933qeZmYGfQ+zuA77j73Ar8XUJSFV/fBZ5Q1ZXAmTj7n7LnWkTmATcBq1X1dCATuI7UPNe/AC4ftyzcub0CWO7+uxH4UTwSZEEncucDe1V1n6r2A/cD1yQ4TXGhqg2qusV93YlzE5qHs7/3uKvdA7wnMSmMDxGZD1wF3O2+F+AS4EF3lVTc5xLgbcBPAVS1X1XbSPFzjTOBZb6IZAEFQAMpeK5V9Y9Ay7jF4c7tNcAv1fESUCoic7xOkwWdyM0DDge9r3OXpTQRqQbOBjYCVaraAE5gAioTl7K4+E/gs8Cw+34m0Kaqg+77VDznS4BjwM/dYsW7RaSQFD7XqloP3Akcwgk27cBmUv9cB4Q7t77c4yzoRE5CLEvp9uYiUgT8FvhnVe1IdHriSUSuBppUdXPw4hCrpto5zwLOAX6kqmcD3aRQUVoobh3GNcBiYC5QiFO0NF6qnevJ+HK9W9CJXB2wIOj9fOBIgtISdyKSjRNw7lXVh9zFjYHstvt/U6LSFwdvAd4tIgdwik4vwcn5lLpFMJCa57wOqFPVje77B3GCUCqf68uA/ap6TFUHgIeAN5P65zog3Ln15R5nQSdyLwPL3RYuOTgVj48mOE1x4dZl/BTYqarfDvroUeB69/X1wCN+py1eVPVWVZ2vqtU45/YZVf1rYANwrbtaSu0zgKoeBQ6LyAp30aXADlL4XOMUq60RkQL3Wg/sc0qf6yDhzu2jwN+4rdjWAO2BYjgv2YgEURCRK3GefjOBn6nq1xKcpLgQkYuA54DtjNZvfA6nXucBYCHOD/d9qjq+knLaE5Ea4F9U9WoRWYKT8ykHtgIfVtW+RKbPayJyFk7jiRxgH3ADzgNpyp5rEfkK8AGclppbgY/h1F+k1LkWkfuAGpwpDBqBLwG/I8S5dQPwD3Bau/UAN6jqJs/TZEHHGGOMX6x4zRhjjG8s6BhjjPGNBR1jjDG+yZp8ldRUWlqqy5YtS3QykkJ3dzeFhYWJTkbC2XEYZcdilB2LUZs3b25W1VmxbCPpg46I/AwIdNw7PcTngjN21JU4LS4+GhjCZSJVVVVs2uR5w4xpqba2lpqamkQnI+HsOIyyYzHKjsUoETkY6zamQ/HaLzh5wLpgvgxSZ4wxJnZJH3TCDFgXzJdB6owxxsQu6YvXIhBukLqTetKKyI04uSHAyTYb6OrqsmOBHYdgdixG2bHwVioEnYgHqVPVu4C7AERErZzWYWXWDjsOo+xYjLJj4a2kL16LQFoNxGmMMdNZKgQdXwapM8YYE7ukL14LHrBOROpwBqzLBlDV/wYex2kuvRd3kLrEpNQYY8xkkj7oqOoHJ/lcgU/4lBxjjDExSIXiNWOMMdOEBR1jjDG+saBjjDHGN57V6YjIPGBR8Dbd0QSMMcYYwKOgIyJ34Ez9ugMYchcrYEHHGGPMCK9yOu8BVkz3+cSNMcbEl1d1Ovtw+84YY4wx4XiV0+kBtonIemAkt6OqN3m0fWOMMSnAq6DzqPvPGGOMCcuToKOq94hIDnCKu2i3qg54sW1jjDGpw6vWazXAPcABnKkGFojI9dZk2hhjTDCvitf+A3iHqu4GEJFTgPuAcz3avjHGmBTgVeu17EDAAVDVN7DWbMYYY8bxKqezSUR+CvzKff/XwGaPtm2MMSZFeBV0/gFneoGbcOp0/gj8l0fbNsYYkyK8ar3WB3zb/WeMMcaEFFPQEZEHVPX9IrIdZ6y1MVT1jFi2b0BEAHDmqjPGmOkt1pzOze7/V8eakESxm7oxxvgnptZrqtrgvvxHVT0Y/A/4x9iTZ4wxJpV41WR6bYhlV3ixYRG5XER2i8heEbklxOcfFZFjIrLN/fcxL77XGGOM92Kt0/kHnBzNEhF5NeijYuCFWLbtbj8T+CFOUKsDXhaRR1V1x7hVf6Oq/xTr9yULK/IzxqSqWOt0/i/wB+AbQHAupFNVW2LcNsD5wF5V3QcgIvcD1+BMFmeSgAVIY0w0Ygo6qtoOtAMfBBCRSiAPKBKRIlU9FGP65gGHg97XAReEWO+vRORtwBvAp1T1cIh1EJEbgRsD72trawn1Oln4lb6urq6Yt5+Mxy9aXhyHVGHHYpQdC2+JF0+oIvIunD46c4EmYBGwU1VPi3G77wPeqaofc99/BDhfVT8ZtM5MoEtV+0Tk48D7VfWSCLatqpqUT+rBafIjfbW1tdTU1Ezpb5Px+E1VLMch1dixGGXHYpSIbFbV1bFsw6uGBLcDa4A3VHUxcCke1Ong5GwWBL2fDxwJXkFVjwdNk/0TUniQUREZuckbY8x05FXQGVDV40CGiGSo6gbgLA+2+zKwXEQWu/P1XMe4yeJEZE7Q23cDOz34XjOBiYKfBUZjzES8GnutTUSKcMZcu1dEmoDBWDeqqoMi8k/Ak0Am8DNVfV1EvgpsUtVHgZtE5N3u97UAH431e40xxsSHV0HnGqAX+BTOCNMzgK96sWFVfRx4fNyyLwa9vhW41YvvMqkvleqgjJmOvBrwszvo7T1ebDOd2I0wfuzYGpNcYu0c2snYgT7FfS+AqmpJLNuPp1MAamrYEFiQwNYp49MQ/D7cay+d1dYGpaURrz9RmpLheAaLJn3RHodUZsdilB0Lb8XaT6fYq4QY47faZ58FoObiixOcEmPShyf9dABE5CJguar+XEQqgGJV3e/JxuMgVD+dRBXFjP/ecGmKV/qi7YcwUZqSrTgr0mML1h8jmB2LUXYsRiVNPx0R+RLwb4xW6OcAv/Zi22b6subT0bNjZlKdV/103ovTR6YbQFWP4Az6aYwxxozwKuj0q1NGoQAiUujRdk2SsCdwY4wXvAo6D4jIj4FSEfnfwNM4Q9IY47vpFiCnW3qNiYVX/XTuFJG1QAewAviiqq7zYtvGGGNSR8xBx51o7UlVvQywQGOMMSasmIvXVHUI6BGRGR6kx0TBimWMMdONV2Ov9QLbRWQdbgs2AFW9yaPtG2OMSQFeBZ3H3H8mQsnWidIYY/zgVZ3OWlX9sAfpMSbl2QOHSWcxBx1VHRKRWSKSo6r9XiTKGL9ZIDDGH14Vrx0AXhCRRxlbp/Ntj7ZvjDEmBXgVdI64/zKw4W+MMcaE4VXn0K8AiEix81a7vNiuMcaY1OLVKNOni8hW4DXgdRHZLCKnebFtY4wxqcOrsdfuAj6tqotUdRHwGTwae01ELheR3SKyV0RuCfF5roj8xv18o4hUe/G9xhhjvOdV0ClU1ZGZgFW1Foh5pGm3OfYPgSuAU4EPisip41b7O6BVVZcB3wHuiPV7jffiPXqCjc5gzPTgVUOCfSLyBeBX7vsPA17MGno+sFdV9wGIyP3ANcCOoHWuAb7svn4Q+IEEpgWdRPBNKtzreJvoeyNNXzKkd6pp91qypy9cGhKZDmP85FVO52+BWcBD7r8K4AYPtjsPOBz0vs5dFnIdVR0E2oGZoTYmIjeKyCYR2eRB2owxxkTJq9ZrrUA8xlkL9bg3PgcTyTrOQtW7cOqfWLFihe7evTvyhAR1HhzfkTDcZxOt57WppG+yz+JpKscs3mmtra2lpqZmwrR6kfbpINyxSEd2LEZ5kQP3qvXaOhEpDXpfJiJPerDpOmBB0Pv5OP2BQq4jIlnADKDFg+82xjOqOm0CjjHx5FXxWoWqtgXeuDmfSg+2+zKwXEQWi0gOcB3w6Lh1HgWud19fCzwTSX2OMcYY/3kVdIZFZGHgjYgsIkwRVzTcOpp/Ap4EdgIPqOrrIvJVEXm3u9pPgZkishf4NHBSs+pkYU+7xph051Xrtc8Dz4vIs+77twE3erFhVX0ceHzcsi8Gve4F3ufFdxljjIkvrxoSPCEi5wBr3EWfUtVmL7ZtImM5KGPMdOBVTgfgzTg5nIDfe7htY8awIGvM9ORV67VvAjfjdNrcAdwsIt/wYtvGGGNSh1c5nSuBs1R1GEBE7gG2Ard6tH1jjDEpwKvWawClQa9neLhdY4wxKcKrnM43gK0isgFnhIC3AZ/zaNsmBaVDnUw67KMx0fKq9dp9IlILnIcTdP5NVY96se1kYTeQ9GXn3hjveNWQYL2qNqjqo6r6iKoeFZH1XmzbGGNM6ogppyMieUABUCEiZYwOvlkCzI0xbcYYY1JMrMVrfw/8M06A2RK0vANn8jVjjDFmRExBR1W/C3xXRD6pqt/3KE0mQazuwhgTb161XmsXkb8Zv1BVf+nR9o1Jeha0jZmcV0HnvKDXecClOMVtFnRMyrHgYszUedVk+pPB70VkBvArL7ZtopeqN8VU3S9j0omXIxIE6wFOidO2jTHGTFOe5HRE5H8YnbQtE1gFPODFto0xxqQOr+p07gx6PYjTX+eDHm3bGGNMivCqTudZETkL+BDwfmA/8Fsvtj2dWR1EZKbzcZrOaTcmEWIdkeAU4DqcXM1x4DeAqOrbY02YiJS726sGDgDvV9XWEOsNAdvdt4dU9d2xfrcxxpj4iLUhwS6c5tHvUtWL3A6iQ7EnC4BbgPWquhxY774P5YSqnuX+s4BjjDFJLNag81fAUWCDiPxERC5ldPy1WF0D3OO+vgd4j0fbjZmqWrGKT+xYG5NaYgo6qvqwqn4AWAnUAp8CqkTkRyLyjhjTVqWqDe73NACVYdbLE5FNIvKSiCRNYEom0+nGPZ3SaoyJnlcNCbqBe4F73bqY9+EUhz010d+JyNPA7BAffT6Kr1+oqkdEZAnwjIhsV9W/hPm+G4EbAWbNmkVtbW0UXzO54O15vW2vBdK3YcMGurq6kj69frDjMMqOxSg7Ft6SZH2qFJHdQI2qNojIHKBWVVdM8je/AH6vqg9Otv0VK1bo7t27vUor4DylB79ORqHSV1tbS01NTYJSlDzsOIyyYzHKjsUoEdmsqqtj2Ua8RiTwwqPA9e7r64FHxq8gImUikuu+rgDeAuzwLYXGGGOiksxB55vAWhHZA6x13yMiq0XkbnedVcAmEXkF2AB8U1Ut6BhjTJLyakQCz6nqcZzm2OOXbwI+5r5+EXiTz0mbULIWqwUke/qMMaktmXM6xhhjUkzS5nSmE8s9GGNMZCynY4wxxjcWdIwxxvjGgo4xxhjfWNAxxhjjm6QdkSDeRKQT8GZIgumvAmhOdCKSgB2HUXYsRtmxGLVCVYtj2UA6t17bHetwDqlCRDbZsbDjEMyOxSg7FqNEZFOs27DiNWOMMb6xoGOMMcY36Rx07kp0ApKIHQuHHYdRdixG2bEYFfOxSNuGBMYYY/yXzjkdY4wxPku7oCMil4vIbhHZKyK3JDo9fhKRBSKyQUR2isjrInKzu7xcRNaJyB73/7JEp9UvIpIpIltF5Pfu+8UistE9Fr8RkZxEp9EPIlIqIg+KyC73+rgwXa8LEfmU+/t4TUTuE5G8dLkuRORnItIkIq8FLQt5HYjje+699FUROSeS70iroCMimcAPgSuAU4EPisipiU2VrwaBz6jqKmAN8Al3/28B1qvqcmC9+z5d3AzsDHp/B/Ad91i0An+XkFT577vAE6q6EjgT55ik3XUhIvOAm4DVqno6kAlcR/pcF78ALh+3LNx1cAWw3P13I/CjSL4grYIOcD6wV1X3qWo/cD9wTYLT5BtVbVDVLe7rTpwbyzycY3CPu9o9wHsSk0J/ich84Crgbve9AJcAgenO0+JYiEgJ8DbgpwCq2q+qbaTpdYHTfzFfRLKAAqCBNLkuVPWPQMu4xeGug2uAX6rjJaBUROZM9h3pFnTmAYeD3te5y9KOiFQDZwMbgSpVbQAnMAGViUuZr/4T+Cww7L6fCbSp6qD7Pl2ujyXAMeDnblHj3SJSSBpeF6paD9wJHMIJNu3AZtLzuggIdx1M6X6abkFHQixLu+Z7IlIE/Bb4Z1XtSHR6EkFErgaaVHVz8OIQq6bD9ZEFnAP8SFXPBrpJg6K0UNz6imuAxcBcoBCnGGm8dLguJjOl30u6BZ06YEHQ+/nAkQSlJSFEJBsn4Nyrqg+5ixsD2WL3/6ZEpc9HbwHeLSIHcIpZL8HJ+ZS6xSqQPtdHHVCnqhvd9w/iBKF0vC4uA/ar6jFVHQAeAt5Mel4XAeGugyndT9Mt6LwMLHdbouTgVBA+muA0+cats/gpsFNVvx300aPA9e7r64FH/E6b31T1VlWdr6rVONfBM6r618AG4Fp3tXQ5FkeBwyKywl10KbCDNLwucIrV1ohIgft7CRyLtLsugoS7Dh4F/sZtxbYGaA8Uw00k7TqHisiVOE+0mcDPVPVrCU6Sb0TkIuA5YDuj9Rifw6nXeQBYiPOje5+qjq9MTFkiUgP8i6peLSJLcHI+5cBW4MOq2pfI9PlBRM7CaVCRA+wDbsB5KE2760JEvgJ8AKe151bgYzh1FSl/XYjIfUANzsjajcCXgN8R4jpwg/IPcFq79QA3qOqkA4KmXdAxxhiTOOlWvGaMMSaBLOgYY4zxjQUdY4wxvrGgY4wxxjcWdIwxxvjGgo4xPhORrkSnwZhEsaBjjDHGNxZ0jEkCIvIud76WrSLytIhUuctnuXOYbBGRH4vIQRGpSHR6jZkqCzrGJIfngTXugJv344x+DU6P8GdU9RzgYZxe4cZMW1mTr2KM8cF84DfugIo5wH53+UXAewFU9QkRaU1Q+ozxhOV0jEkO3wd+oKpvAv4eyHOXhxo+3phpy4KOMclhBlDvvr4+aPnzwPsBROQdQJnP6TLGUzbgpzE+E5Fhxs478m3gL8B3cALPS8B5qlojIpXAfTjB5lmc0Y8Xp+IIxyY9WNAxJomJSC4wpKqDInIhzuyeZyU6XcZMlTUkMCa5LQQeEJEMoB/43wlOjzExsZyOMcYY31hDAmOMMb6xoGOMMcY3FnSMMcb4xoKOMcYY31jQMcYY4xsLOsYYY3zz/wGh1xazr/3i0gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1e274103a58>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "x = 2*np.sin(np.arange(100)/2.0) # periodic signal\n",
    "x += np.random.randn(len(x)) # corrupted with noise\n",
    "\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot(x, '-s')\n",
    "plt.ylabel('x', fontsize=20)\n",
    "plt.grid(True)\n",
    "plt.xlabel('Time')\n",
    "\n",
    "plt.subplot(2,1,2)\n",
    "c = plt.acorr(x, usevlines=True, normed=True, maxlags=50, lw=2)\n",
    "plt.grid(True)\n",
    "plt.axhline(0, color='black', lw=2)\n",
    "plt.axhline(1/np.exp(1), color='red')\n",
    "plt.ylabel('Autocorrelation')\n",
    "plt.xlim(xmin=0,xmax=100)\n",
    "plt.xlabel('Lag')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>图5.17:显示使用正弦函数生成的100个随机数的图，以及该系列的自相关。</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "自相关函数也用于计算相关长度。相关长度是指从这一个点到另一个点的距离，与该点相关的物理属性没有进一步的相关性。在数学上，相关长度是自相关等于$e^{-1}$的lag，这由图中的水平红线表示。让我们绘制一个具有较高周期性的图来计算相关长度。结果如图5.18所示。从图形上看，相关性大约为9。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEKCAYAAADXdbjqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xl8XGd18PHf0TLaN0uybMv7mjgLWZwFAkFOYiCQJk1fCoQuQKHpAg3Qvm2T4IRCgYaSQtnblJ2XJqRhSRqCg5NYhASyyInjLXa8xI53ydpH20ij8/5x78xcyTPSSDOjWXS+n48/1oyuZp57584999nOI6qKMcYYk4i8dBfAGGNM9rNgYowxJmEWTIwxxiTMgokxxpiEWTAxxhiTMAsmxhhjEmbBxBhjTMIsmBhjjEmYBRNjjDEJK0h3AZKtrq5Oly5dOu2/7+vro6ysLHkFygKzcZ9hdu73bNxnmJ37PdV93rp162lVrZ/u++VcMFm6dCktLS3T/vvm5maampqSV6AsMBv3GWbnfs/GfYbZud9T3WcROZzI++VcMDGxrfvMZk77A2c8X+mD7U0zXx5jTO6wYDKLRAskAD3RnzbGjBPrhqyu3EfLxg1pKFHmsA54Y4yJU6wbsljPzyYWTIwxxiTMmrmMMQmz5h9jNRNjTMKs+cdYMJlF6sp9UZ+vjP60McbEzZq5ZpGWjRt46UgXN3z9aQB+8ldv4OIlNTQ3N6e3YMYkSaqHv9eV+2I25812FkxmmRPdA+Gf2/1DaSyJMfGLt08m1cPfWzZuQFVZvfGXDAeV//3IGzlvYVVyXjzLWTPXLHO8azD8c0eftWeb7JBJfTJ9gSDDQQWgZ3B4xt8/U1nNZJY50T1AYb4wHFTaLZiYJMmm5p9ER551er43PQMWTEIsmMwyx7sHWVhTSmvPIO020sYkScvGDU5fnCrbj3XzsatX89FrVqW7WFElWsvp9gQQq5lEWDCZZU50DTC/qpjgqNLeZ30mJjlUlYOtfm68qJH+QJCXjnalu0gp09nvrZmMpLEkmcWCySxzvGuQK1bW0R8IWp+JSZrW3iF6h0ZYObec/kCQLXtaUVVEJLzNTExsjNXclszh7539VjOJxoJJhkrFF28kOEpr7yALqovpHghwzNMZb0wiDrT6AVhRX46I8MDWoxztHGDRnNLwNok0L8XbJ9OycQN/+cOtbNp1EoDH/vZKVs6tSOrw9263ZiJifSZeFkwy1ERfvKW3/uKM5+MJMqd6hxhVmF9VwqmeQXYc6552+Sx9hvHa3xYJJpXFhQC8dLRrTDBJRMvGDXT1B7jg05sBeNOqOn74wcuibnuiZ5C5FUW09g6x41g3K+dWJKUMIaGaSUNFMT2D1swVYkODc0Q8d3cnupw5JvOri6ktL6KjL4CqJvX9LH3G7HSg1U95UQENlUWsmVeBryCPl44kt9/kcHs/AOVFBRxq74u53alupym3uDCPHUd7zvh9rBFm8Y486+wPUF5UQG25z2omHhlfMxGRRcAPgHnAKHCPqn45vaXKTse7nWatxuoSast8DAf1jDsrq3GY6TjQ1seK+jJEBF+BcM6CSl46Mv2abzSHO5xgcsXKWjbvPkVgZBRfwdj74VBT7sKaEtbOr2RnlNp3y8YN3PbTHTy8/ThDw6N84Iql3Pb2s+MuR3f/MFUlhVQWF1qfiUfGBxNgBPg7VX1BRCqArSKyWVV3p7tg2SZcM6kqpta9Cxs/C95qHGY69rf6ecOKWmDsDUm0JtloQttNdNPymlsbedOqeh7ddYojnf2sqC8fs02b32nKnVdVzHmNVTyw9Sijo2fWvg+2+VndUEGewDOvdsS3k67O/gA1ZYVUlhRw6HT/lP42l2V8M5eqnlDVF9yfe4GXgcb0lio7negepKKogIriQuaUFQE2C94kzj80wsmeQVbMdS7sidx4TPS3h9r7mVtRxNnzKwE4HKWp66Rb+55XWcy5jVX0BYIcPH3mdq+e7mNZXRmXL69l57Fu/EPx93109g9TU+qzmsk42VAzCRORpcCFwLPpLUnyxWpeSqbjXQPMry4GoLbMqZmc9gcoTum7mlznHck1mSf/fj1XfmHLtN7ntfZ+ltaWsbTW6dSPVisIB5OqYhprSgDYeaybas82vYPDtPYOsby+jPMbq/nqE/tpOdRB05q5cZWje2CYhTUlVJYUWp+JR9YEExEpB34CfExVe8b97mbgZoCGhoaEhgH6/f60ZNGdKJBU+iZPVFfpY9Jy7z06QKVPaG5upmNwFIBnXtzBupqhuPbZu02sMsVTjkyRrs86nVKxz08fcy6oHYd203x6z4TbPvHU7yZ9vVjle+VEP+fW5bP9+d9SUgBPb3+F5SOHx2zzm0NOWQ7ufIHSAvDlwSPP7OT6RcPh1z3UHQSg/+Qh+kbyyRe4v3kbnIivA761u49lpUN09At9gSCPP7GF/DyZ/A9n2Eyf31kRTESkECeQ/EhVfzr+96p6D3APwLp167SpqWna79Xc3Ewifz9tm2K3LW//9DvOeO5Am5+r/+3X/Os7z+dd6xbF9Rb+32zmDcsbaGo6n6GRIH/bvInaBUspzz8W2ecJyuE9Lvet7uHtX/kNeQLnLazmwQ9fEVcZMknaPus0SsU+P//oHgryDvKH1zZRmJ834Tm06pzXwe8mbliIVr6BQJCuTZu4bO1y1q9fxcqdTzFc7KOp6dIx2/3ukZfx7T/EdRuaEBHO2fM0XZJHeflQ+HUf3HYMfreN69dfxuqGCi7Y91tOBJWmpsnP4eCo0v/oI5y7cik1ZT5+tn83F112BTVlmZeDbKbP74wPJuJMof028LKqfjHd5ckUy2rLKC8qYMfR7riCyeBwkPa+APOrnKp/UUE+FUUFTrLHysh28U4O2/paJwCvX1HL3pP+BPbEzKRbnuijJ8rFfjqj9cY3za76xC8n/ZvuaTYLveaO5FpSV+b8X1sadZ7Uie5B5lUWIyJjyvfcIcJBrtSXjwgsnlMac7BArOPROziMKlS5fSbgzILPxGAy0zI+mABXAH8C7BCRbe5zt6vqI2ksU9rl5QnnNlayPc6Jh6G25PlVkR6S2nLfGcGkZeMG9p3qZcOXnuSdFy/kga1H+fb71nH12Q1jXu+Fw53MrSjismW1PL2/ncHhIMWF+YnvmEmpWM2l0+mvm+rf1JX7wsFkTpkv6uCPWHM9Qp3tS9xJkMvqyvjlzpMMB0ed2pDrZPcg89xzPFb5+gNBFs0pobgwf8qjF0MTFmtKC6kIBZOBERtSTxYEE1V9Csi8BskMcP7Car7320NRx9uPd9xdFKuxuiT8nPOFPjPZ4yunnJrGH122mF/uOMGWva1nBJOthzu5eElN+PWOdw2wPI4OWDM7HLrrHdxy74tsO9LFk/+wPvz8f/76AABP/sN6yosKUFXO/9SvuP51C/jsjefFfL1wzcTtfF9SW0ZwVDnaOcAyt7YCcKJngIsW10xavmV10ztXQ0kea0p9lPqcm6eewWEbUk8WBJPZYjrrQZzXWEVgZJRXTvVybuPEq72dcPNwzfcEk9ryIo50nDki5pVTveQJnD2/kitW1rFlT9uYpH2tvYO81tHPn1y+JDxi5njXoAUTM0ZNaeGYDLvgNHPl5wll7oVYRDh3QVV4cmGsO/ziwjyqSgqpLnW+D5ERXX3hYKKqnOoeCtdMJrLcE4CmotutmVSVFlLi1sRtRJfDgkmGaNm4gbt+uYdv/eYgL965IVyFnsjrFjoDHncc644ZTMZ/Odff3Qw4QeqasxvYdqSL8afBvtZeFs8ppbgwn/VnzeVXu0+xv9XPqgYnx9ELh500GRctqWFuhTNf5ViXTd4yY9WU+egdHGEkOEqB2xTVPeDMHvdmEz5vYRXf++0hhoOjMe/kB4dHWd0QybG11A0G3rQqHX0BAsFR5lXGEUzqpxdMvDWTUGvAZHNN4pmQmQssmGSQ5r2tXLJ0TlyBBGDRnBKqSgrZfrSbmy6Nvs1E1e/acqfdelTH9nW8csoJHN5AtOFLT4Z/X1KYjy8/j3MbK8kTIU/gWOcAxnjVuLWIroFh6sqdm46ewREqi8deds5ZUElgZJT9rRMP5FjsSRpZW+ajvKggnK8LnM53GNsvGMvyaTZzdXn6TELDgeNd0yTXm7wsmBClar1p5u8kjnUNsOdkL7e//ay4/0ZEOH9hFTuOTS+h3pyyImeoo+fGKjAyyqHTfbz1nAY27z4V9e8Ghp1x+ms2bqKu3EdDZTFHuyyYZINY84Oms7zuZE2z1aXOTVFXfyAcTEI1E69QrXqyLNZLayO1CRFhSW0pr3pmt5/qCU1YLJmwfBCpmUy1ebmrP4AIVBQXIkCe2JomIRZMyIx8VM17WwG46qz4ZuFC7BxI8QbB0BemJxDJXfTq6T5GRnVMk8JETvsDrFtSYzWTLPGVq8q4/Io3cdYdmwD4gwsb+eK7L5jWa7Vs3MCdD+7k5y8eY/s/vfWM34dqJt7FpLoHhqkcF0yW1ZZR5stn1yTB5Gtb9vO1LfsRwJttK3Tulxc5NexQzcT7Hfj0/9vMd3YGuHhJDbuOd4ebwrzb/PCZw9zx8508+rErWTMv+vnf6SZ5DNVKKoqdWfATBa7ZwoLJNCRrGGC017nmi0/G/TqJBsE57tj4Xk8weeVULwCrprAGRGNNCVsPd8a9fTxsqGXqhIaJA+ybpGlpMie6B8Nzl8YLBxPPEOBeNxWJV16ecM6CKnYePzNdfDSxFk3wDwXJz5NwLcjrgVecMoTO0+W3OzMLvOfT286Zxycf3MkvdpyIGUy6BobD+wVQWVJAz+AILRs3oKqcdccmhkZG49qPWLL13LdgMg3JqsmkskYUT7bWWjfZo7dmss8dyTWVDsrG6hJ+sf0EwVFNWlqJTKgt5qrjbpPk2fMr2d/qZ3RUyZvm53aiO5LvbbxIM9fYmsn4Zi6Acxorue+5I9SWuXOfpmluRVHUczCe+TXXfvlJRhW+8vg+vvL4vvDz3ot4V39gTPkriyP5uU77AwkHkvFliuf5TJHxWYNNatSV+8LNXGNrJn6W1pZNaQJiY00JI6NKa68tA5wNjrnB5M2r6xkYDoYfT8fJ7sGYHd6hWeGhEVCqGjOYnLugioHhIO+9bDEA//2h6KsoTiaeYcGxxHMR7+wPUFM6Lpi4fSah4zh+gEHIdPqlsonVTCYR73oM2eLQXZE8X8NB5y5qTDBp7WWlm0o83nbg0MTFY50DMZs8TOYIjXq6clUd//HrA+xr7R2zvG68zSxDI0FO+wPMq4z+mZf58inMl3CfycBwkJFRDachifZ+X31iPwDv/db0EoPHM5IrEV39w6z2NAFXlRSGBwGE+g3vu/n1rF3gpJXYsreVD3z3ef77Q5fxhpV1KS1bulnNhNy/Y4ilMN+ZCBZq5hoaCXK4vT/c+d6ycQOH7noHh+56x4RLnYbawBO5wzUz53jXAPUVRZyzwBlFFcp4EBJvM0trj5M9IdYFXESoLvXR5dZMQqlUxtdMktl80xDHHJNEdPUPU+WtmZQUeGomzjDlRk+f0EWLnNn4ye5TzERWMyEyouP9332OQyfaab79WiDzayWJjCDx3g0+/trYfQ2NmvHeiU7U8dcfcMbZH7URXVnhWNcAC6pLqCotZG5FEftOxd8J7z1PQkEhVp8JjJ0FHyuYTNX40Vxe3336EN99+lBKOqsDI6P4h0bGdsB7+kyOdQ5QUVwwZv+qSgtZ3VAeToyayyyYeNSXF/HSUKzTdHKx7t5jNRvE+lLEW1PyflmmGvjiCULxBqpSXwE1pYVJrZlMJ72Mic/xroHwaKVVDeXsb+2d1uuEgsNETUvVpb5wM1docl9lSfyXHW+zbDSxzvvx504y5teE9ndMn0lJIX2BICPBUY52DozJfRdy8ZIafrH9RNwDHbL13Ldg4lFXUURPQOP60J++9SquuOsJADa+42w+9KblMbeNdVH2BpJE76TSPc69saYkPEooGVo2biAwMsrqjU5a8y+9+3XceOHCpL3+bKWqHO8aCq8quGpuBfe3HEloRNe8CfrJakojfQrJqplMx1euKpt0bY/JLuKh5rqqMTUT5xLaOzjCsa6BM4Y9A1y0uIZ7nzvCgbZISqIJy3rThbz3v57lmrMbeOzlUzz59+tZXFs66d+lmwUTj/ryIoLqnPSTrU/gHTu/96RzZ5fI0ruJBgJvIJqoAzVVGqtLONh25lrbiejyJAkMtc+bxPQNOx3hC9w76FUN5fQHghzvHmBhzdQvWBXFBZQXxb6M1JT6eKHfydCQzmASD+936O5H9/LNXx9g68ZrwsklvennQ0ITMHsGhznWOcBly+ac8boXL3H6TVoOd8YVTLbsacWXn8efvH4Jj718igNtfgsm2abeTVrY5h+ipsw34Z1KaCx8mS8/PNEvU8aBp2NiU2N1Kb/Zd3pMduFEdXiDSa8Fk2Rod5drbnT7OUKTU/e1+sPBZCq13MlGT4U64EPDgoEzRnNlYrPOhrUNfG3Lfp7Y08ofXOTUiLs8SR5DQvtytHOA3qGRMZ3vMPbG7raf7uC2n+4AJm6J2LK3jcuWz+E8N83MgTY/66eQGSNdLJh4hINJ7xCrGypo2biBbzTv51837SU/T9j3mWvDTQE/f/EYAJcum8MzBzsYHZ1+X0u2835hlt0WWbMsVp/QRH1F3i+Yd/EkCybJ0THoHPlwzcQdBr7vVC/r3aavlo0buP1nO3ho23G23bmBlROsoDhRExc4d/HDQaUvEAx3VI9Pp5KJs7rPa6yiobKIzbtPeYLJmTWr0L68fMKZvT++djfVCYhHOvrZ3+rnvZcuZk6Zj+rSQg64Nf5MnxlvwcQjlIahzXPhOuWOyQ+O6pjmr1DN5PUratmyty28eE+2iOfuM947w3j6hOJ5fvzrdPZFvrytPTYhMhnaByLBxHtx+twje/jcI3sA53MvKyrg8uVzKMjPi3muCLBgkpqJN6VK98AwFUUFScuSECprKmo1l37uMU77A/xy58kzOvm9TeChwQS73WASrQM+HuMDxacf3s2nH95NQZ5woM0ZbZfpM+NTEkxEpEBV48vLnEFCNZPT/kgwOem5iLX3BSKzevsC5OcJlyx12kj3nJzeiJh08d7JNDc3T9o5mQ6hVSDXzKsYE+BNxFTvVtsHFV9BHrVlsW8mTvsDnPYHeP8blgKRc+XBbcf46H3beOSWN7FybjmrN/5y0hnn3pQqPVGSPCYqVXfkE12gz/3ko4BzjB/6yBsB2O3mFRvfzJXo+42MatL7IlNlSpMWReQeEZnw7BGRZcBTCZUqTSqLCyjIG1szOdkzhM9d2KfdE2Ta+5y0CmvmVSASSZA4kTuvWzvh5D/jDPVceusvWPeZzXS4NZOz5lWE04ubsaZ6t9o+MMqCquK4+rWuGDdjO7QY2/ajXeHPY7I+E29KlZ7B5AeTdDrtD4T350Cbn+JCJ0gn/32Gwv1NmWyqNZMPAZeLyLtUdc/4X4rIO4H/AiqTUbiZJiJU+eSMZq418yrYcax7TAK6zr4Ac8p8lPoKWDynlL0neydtOvr0w7upK/dNOnbeOF/Uzv4AFcUFLKguoS8QpG9ohLIJRg6ZyXUMKgtq47t7fsuXxmawXlJbSmVxAduPdbPC7WuJp88EnGDi5OXKrc+vzJdPnsBwUFk0pzRpg0/GO9iWWHbnmTDVdCqfBdYCLSLygdCTIuITkW8APwaCwI3JKyKIyNtEZK+I7BeRW5P52uNVFQltbg0kOKq0+YdYO9+Jjd5g0tEXCLcHr2moYM/JHn7211cAzryTWDKlfTMbdLgBO7Q0sHXCJ659QMOd7/Hwnq/OYmzV7DjaHZ5TNFmfSWhYbVf/cMwkj9lMRMK1k2j9JclqiTiQBU1dUwomqnoH8FagF/iWiPxQRNYBzwF/CfwWuEBVH0pWAUUkH/g6cC1OILtJRNYm6/XHqyqK1ExO+4cIjipnz3eGT3qbuTr6nWVvwWnTP9Tez6adJ4GpLXCVC6b6xYj33s3J0OpjboVzwbJO+MQMB0fpGppaMBnvvIVV7DnZw2vucrmT9pmURGomPQMjZwwLzgWhfYo2T8eb3+73L1jAgqpiDt31jin19RTkCQfb/BnfRD7lRI+q+jjwOuAx4L3As8A5wGeAN6vq0aSWEC4F9qvqQVUNAPcBNyT5PcKqfBLugA8tItRYU0p1aSHt/rHNXOGaybwKgqPKd59+lWV1ZSyvn9760tnK+4WZyGdvPJdDd72DV+96BxcvqeGSpTUTbh+umVRazSQZTnYPokTmmEznInR+YxXDQWXL3lbKiwqomCQ4FOTnUVFckHU1k6kcm9CIrmiz371W1JdzvHswnMsunverK/exuLaUA21+fvShywH43I3n8dGrVyECz95+dUYMC4bpj+byA21EbjK7gSdVNfGVYc7UCBzxPD4KTG+xgzhUFgntfQFGgqPhkVzzKovdRXuci9noqNLZ71zovKNpjrvBJ9MTRKbLzmPOiJfgqPLyiR7etW4Rzx+KnQCvsy/A2fMrrZlrAlMZGhtKPR+qmUwnt9v5i5xO+Bde6wrPUZlMTamP1t5BBoaDWRNMppJRIlQzmWxYcOgm82BbX3jde+/7/eevD/Avv9zDrk+9dUzf4J//oIWDbX08uuskInDN2rn0Do7w5cf38fD2E3zwjcumt5NJNuVgIiKvw+kbWQU8CvwMuBvYJCKfB+5IclCJ1ioyZqqCiNwM3AzQ0NBAc3PztN+smACqwsOPNdNyMgjAgV1bKRgZ4sDRAZqbm/EHlFGFjuOvTasPJJHypYLf709amWIl1MsXeGbPUZqb2znuH6U/ECS/53jM7St90NY7iL/9FC8+20FBHrTs2seKkcNJKSckd7/T5e43FrLzdB53twySL/Ctt0Q6gUP7dssTfWOO8Z98+znAOcZfuWryFTW9x0hVqfBBbwCKggNxHb/84CC7DjuB7NTRQzQ3H4tv55Iokc/67jcWAmcGwVue6BsThD/242187MfbYh7Xjl7nsvjwk89zev6Zl95tewMU5MFzv/3NmI78woEAB9uG+Z9n9rOiKo/dW5/hliecPpR/fng3//zw7vC23vee6fN7SsFERD4MfMH9u9tV9fPu81twAsytQJOI3KSqryWpjEeBRZ7HC4Hj3g1U9R7gHoB169ZpInMmWk4+Bgyx6ryLeUVPULD3INdvWM+jbS+wr9VPU9Ob2d/qhyd+zaUXrOVHe7ZN6fXryn0ZN6cjmfNMtsd4mc898jLfe/oQV7zpSh7ZcQLYxjuvvow7/tgZ3NByqIN3/sfv+N4HLqFpzVz6AyOsvfNRXnfWCtY3rWDec09QXD2HpqYLklJOyNz5NVPV+eJRaHmJoMJlV7yJUt/Yr3XPpui1jp4A4f2veyr23bf3GK37zGZ63c12tgd5/6a+8Haxmlu+e/A5dz2PUdadv5amCxuntoNJkIrPOp7j6jU4HOTO326iqG4xTU2rz/j9L9peor7jNOvXrw8/59SKnGHBR9xgFDrmk733TJ/fU62ZfBV4DbhJVX8XelJV94nI5cC/AR8GXgRqk1TG54FV7vyVY8B7cPpqUqKqyLkjaOsd4mTPIHMrisjLE2rLfTxz0GlmCa3PMGeSMeU2BDjinAWVBIKj7G/1s+t4D76CvPCKjhBpejne5dzBhlKpzClz7gjnVhTZssAxePvyugeGzwgm8Yi33X06s7BrSgvxDzn9BNnSzJUKxYX5LKwpCc9oH887QjQkm0Z/TrUD/kHgQm8gCVHVgKr+DfAHSSlZ5HVHgI/gNKm9DNyvqruS+R5eoWBy2h/gVM8gDe5oldqyIroGhhkJjoa/vOM/eBNbaFW/nce62XW8mzUNFRTmR06/hspi8vMkvFpdKJVK6BjPrSi2zMExtHlGGWbi5LZqb2LEHJtnMlUr6stjzmj3jhDNRlMdGnyjqk64ZJiq/hxIXluE85qPqOpqVV2hqp9N5muPV+nz1Ey6B5nnLgNaW+5D1UlDHaqZZPMHP9OW1ZVR6stn1/Eedh3v4dzGsfNa8/OEeZXFkZrJuNrf3Moi64CPYUzNpD/zgon3pms210wAlteVc/C0P2pi2M4oNZNskpLbBFU9MvlWmam4QCjz5dPWO8SpniHetKoecGomAO19Q+EmmJrSidPUm4j8PGHt/Eo27z5FV/8waxdUnbFNY00Jx9ylf0PrxYTSccytKKJ7YJjB4SDFhflJKdMtT/RFbffOlCys8Wr3D1GQJ4yMakbWTGrKzsyyO1utmFvG4PAoJ3oGzxj91e4Ohc9Ws7vOGUNdRRGH2vvwD42EJ2WFaiEd/gAdfQHKfPkUF+Zn1UUn3c5trKLlsFOxPXfBmRl3GqtLeO7VDiDSZ1JbFmnmAqfGuGhOchYKijaKDLKrnRqc8i6uLeVgWx89g2fOYagt942pvYTM1A1PdZT1P3LBdG4kV7jDgw+0+scEk+HgKL2DIwkHk3TexFowiaK+vIidx7oBIs1c7od8us8JJpOtxGjGGj9W/8Zv/BYYWwtorC7hZM8gI8FROvsD5Enk4lMfnrg4OGkwyfR1H5Kt3T/E2gVVHGzri1ozufO6tXz0vm3cd/PlDL62I6ERPtO5gIbycxUV5CWtVpkJpnMuLa93hu0eaPNz5er68PPja+IhEx3v0Pvf9tMd/O9Lx3nhjg34CqY8Dz1pLJhEUV9RFL6Dbgj3mbjNXH6nmSsV2UFzWTyjgBZUlxAcVVp7h9yszL7wYmThiYtxdMJP9F65NqFUVTndF2BZXSkikQ74aAH1Pfc8Q6Uv9vDteEznAhrqB5jt/SXg3KhWFBec0Qkf6iMcf12J53ivX1PPvc+9RsvhDt6wom7S7VMlfWEsg4XWNYFI7qHqkkLyxOnstJpJaoTWgjjeNeB0RnqOcSioWyf8WL1DIwRGRqmvKKKiqCC8mmGsgBqraS+VQmuazPb+EnASQy6vLz9jeHBHAiNEr1hZR2G+0Ly3LSllnC6rmUQRWnERIs1ceXnCnDJn7feOvgCrGmZX/q2ZEMoZdaxrwMnL5flizSn1UZAnNtdknFBfSF15EVWlhRnXAe+tIe1v9Ydrhrna5DgZ7/Hw1pIrip1L8XT0dyMUAAAgAElEQVRGiJYVFXDZslq27Gnl9rfHzlieahZMogjVTCqLCyjxRdp4a8uKws1cc7J4CF+mCk1cPNY1QGd/gGV1kZQUeXlCXXnRjMw1yaaReKFM1rXlRVSVZF4wyfSlZmdarP3udQdOTKdmEitAJdqkOVUWTKKod2sm49Nr15b7ONY1wMBw0Jq5UqDUV0BNaSHHOgfo6Bvm4iXOMfZ+Wf5n61H+Z6uTmDrZd7fZmLEglOG6tsxHZXHmBRMzNaHBClORKU2aFkzGcZLitQDwyqmx1fLXr6jjhdecjnnrgJ+aeEcBNdaUhGsmobu0qd7dTrbiZS4J7Wd9hVMz2d+a+SvymeiqSgopyM/ebmwLJuNMNPegtszH4LCTbM1qJlMTbw1iQVUJ2492ExzVKY25nyxNeDzBZTg4OibFy2Svmwlt/t7UPt5mrlgBtdJO24yVzRMWwYLJlHhrI1YzSY3GmhJ+tfsUMLUv10S1F2/z1URDg9v9gTOaNjO9zf+0f4iqkkJ8BXljgknLxg2c6B7g9f/yBJ+78Tzee9liIPOWPzARFkxmkVrPKC+rmaSGd1ZwPMc4mfNGTvuHJl2GdjIzXZNp7xuKLNJUUsjQyGg45cwpd7DCXM9Q95lm6YbGinU88vMkq/NygQWTKZljNZOU8waTVIyYm6g/xZt9d7pmuiZzujcQvskJTQrsGRimuDCfVnel0NCyx+mQCU2BmSR0PP7poV38T8sRdn7qrYgIl33usWlfUzKlSdOCyRSE7qby8ySncgxlkgXeYOJ+uZLZoe69uIUWD3qtvZ8rv7CF01k4IfJ03xBnz3PynIUmBXYPDDO3sjg8wTM04dNkjuX1ZfQFgrT2DjG3oojOvuFpt3bECtgz3aRpwWScWMvI1pX7wneANaWF4TQfJrlCs+Ah0szl/bJc+a9bOH9hFV9770VJa+Kqq4i/kz7TtPsja2CEayaDTr9Ja88gIlaLzkShOVQH2/oo9eUTCI5m/edkwWScr1xVFjUR3rrPbGb93c3A2BxPmTKqJxeM728495OPAmOP8fL6spiLC0UTT9t8qa+AUl9+eM7G+L/P1Db/wMgo3QPD4eURqjw1E3BSz9SWFWX1cNNctdzNHnzwdCR7cLb3w1owiVOmj+rJBfEc4+V15Tx7sCPq4kJeU52AWFdeFDWYtGzcwMaf7+AnW48xMBzk79+6hg+vXzml106VUJr+UM1qfDA51TNIQxr7S0xs8yuLKSrI49W2PtbOd5op55Rld9O53bKYrLK8voyB4SAnewZj1g6mU2twaiDR+0wOtPaxZl4FC6qKJ50UmMwyTSYy+31czaQ/UjNJ50guE1tenrCsroxXT/eFV26dU5bdn5XVTExWCa0HcbCtj5aNG3jopePccu+LPHLLm1gbZcGteNWVF3G4vT/q70JrT1QUF0waTFo2buBIRz9v+tctAKxpqODRj1857XJNJBRMwkOD3WSB3QNOnqfW3iHOazxzRUuTGZbXl/HyiV46+pzgn+35/qxmYrLKCk9bM8Du4z0U5gsr5yaWxbmuoijq0OCewWFae4dYUV/OCjd1uOrETWxHOpygdMGiag60+QmMjCZUtli8GYMBCvLzKPPl0z0wzEhwlNN+q5lksmV1ZbzW0c8pdwj3nAzoh0tERgcTEfmCiOwRke0i8jMRqU53mUx6za0oosyXH+6Ef/lEDyvnViS8wlxdeRGd/QFGgmMv/KH3WVFfxsq55fQHgpzonjgN/pFOJ5hsWNvAyKiesXZFsrT3hTIGRy5CoVnw7X0BVKHehgVnrGV15QRHle1Hu/C5NwLZLNObuTYDt6nqiIh8HrgN+Md0FCSTR/XkiniO8fjFhXaf6OHKVfVn/M1U1Zf7UHU6ted6LsAH3GatFXPLw/M49rf6uf5rT8Us63suWUx+ntC0pp4vPLqXvSd7OXt+pAkuWbPkT/sD+AryKC+KfI0rSwrpGRwO3+02WM0kY4WabLce7mROmQ+R7J5ukNHBRFV/5Xn4DPDOdJXFhv+mXrzHeHl9GS2HOmnrHaKtd4iz51ck/N6hpqI2/9DYYNLmpyBPWDynNDxRdX+rf8KRZ0c6+1lQXczqhgoK84U9J3vP2CbW307Faf8Q9eVFYy5CoZpJaN2XuVYzyVjL3bkmp/2BMTcb2Sqjg8k4fwb8ON2FMOm3vK6ch146zovucgCJdLyH1Ll38OMv6Afa/CypLaUwP4+6cicz72TNVkc6+llU4/zNivpy9p7sSbh8XuNrNt45TxctruG1jv7w7HfrM8lc1aU+akoL6ewfzvphwZABwUREHgPmRfnVJ1T1QXebTwAjwI9ivMbNwM0ADQ0NCaUR8Pv9sy6zarbt82DbCKrw3ce3AdBxYAfNR6beRODd75N9Tl/JU89vQ49Hvtg7DvUzrywvvF19UZCWV45O+Lr7T3Zzwdx8mpubqZFBXjoc//GNZ7uJajYD3e20dgV5dvseAHa/8AyveLI1ZNtnnSyZut+1viCd/TDi70p6+WZ6n9MeTFT1mol+LyLvA64DrtYYw2hU9R7gHoB169ZptBns8Qrla5pNsm2f6493882XnmJHu7Cgqpjr3rJ+Wq/j3e/ewWFu/c2vqFu4nKY3rwCc9U3aNm/ihkuW0tR0FgCPnH6JJ/a0ArGbpHoCyqVrl9PUtIqXOcDvNu3hwkuvoCq0it6m2GlgvJ9DrL6Viaxetoitba9RWjuf2pMnueaqsccm2z7rZMnE/XY+X+cm5tmTQZ7d5Az2SFZWjZne50wfzfU2nA7361U1+iQAM+uE8hr5h0aS1tZcXlRAcWHemImLRzr6GQ5qeDgywMq55XFd4BfNKQXgLLc/Z++p3ok2j2o62RWqSgrpDwQ51jVg/SUZLteyamR0MAG+BlQAm0Vkm4j8R7oLZNJr3Wc2s/bOR8OPH9/TytJbf8G6z2xO6HVFxE2pEvkiH/AMCw4JzWcpjTGMMzRxMBxM5rnBxNNvkspZ8qFZ8Ada/dZfYmZU2pu5JqKqmZEEyWSMVN7NhfJzjW9euvEbvwVAgFA7a38gGP69ryCP+vIinr71Kr739Kv80//uZlFN6ZjXuePBXdzx4C73fXz84M8u5U+/8xxfvelCPvbjbfzlm5fz9289K+F9CAWTY10DvGFFbcKvZ0y8Mr1mYsyMqSsvoq13KGZgijXvPTAyyrGuAU52D3Kkc4CSwvwJ12A57Q/QcqiDPIH1Z83lkqU1bHaXKk6s/L5wMAFbx8TMrIyumRgzk+orfGw70jXtv996uJMjHf0srCmZdALa84c6WbugkvKiAjasncc/P7ybw+19LKktm/DvQqpKCnnuE1dTVDC2uW3r4c7wz+lcYdHMPhZMjHHVlRfR0Te91RaLC/OcYNI5wGK3v2Qi24508e5LFgHw9Sf2AfDmLzSP2cbbrDZe98AwazZuOmPkT1VJ5CttfSaZLdeyalgwMcZVV17EJMukxHT+wmq2Hu7gSEc/ly2bM+n2A8NB1i2tAaDDTRk/ngIFecLWjRt43ad/FXWb8RejSk8zl43mymy5llXD+kxMVknlSKhQSpXpWLekhh3HuvEPjbDQs/TwxH8zedC5fHltZH5KHLx9JlYzMTPJaiYmq6Tybi4UkAryhJEoVZRYzU515T4uXlITrtWEhgXHasbw5efRUFXEvKrJaw5vPTdacojoxo9Ce+Pnt4TLkWt3wSbzWDAxhrEXYm8giedCvO4zm/ng91vCj//ih1vDfxtaPtj7+oHgKEc6Blh66y8mrVHd8fOd3PHznXHtQ65NgjPZxZq5jCGxC3E8f2sXepPrLJgYk2aJ9Pdk68gfk3usmcuYNPM2o4XSyUcTajIzJhNZzcQYY0zCLJgYkyNSOWzamMlYM5cxJDYbOZ6/jff1EymHDf816SQx1pvKWiLSBhxO4CXqgNNJKk62mI37DLNzv2fjPsPs3O+p7vMSVa2f7pvlXDBJlIi0qOq6dJdjJs3GfYbZud+zcZ9hdu73TO+z9ZkYY4xJmAUTY4wxCbNgcqZ70l2ANJiN+wyzc79n4z7D7NzvGd1n6zMxxhiTMKuZGGOMSZgFE2OMMQmzYOISkbeJyF4R2S8it6a7PKkiIotEZIuIvCwiu0Tko+7zc0Rks4jsc/+vSXdZk01E8kXkRRF52H28TESedff5xyKSc1PFRaRaRB4QkT3uZ/76XP+sReTj7rm9U0TuFZHiXPysReQ7ItIqIjs9z0X9bMXxFff6tl1ELkp2eSyY4FxkgK8D1wJrgZtEZG16S5UyI8DfqerZwOXAh919vRV4XFVXAY+7j3PNR4GXPY8/D3zJ3edO4INpKVVqfRnYpKpnAa/D2f+c/axFpBG4BVinqucC+cB7yM3P+nvA28Y9F+uzvRZY5f67GfhmsgtjwcRxKbBfVQ+qagC4D7ghzWVKCVU9oaovuD/34lxcGnH29/vuZt8Hfj89JUwNEVkIvAP4lvtYgKuAB9xNcnGfK4ErgW8DqGpAVbvI8c8aJ01UiYgUAKXACXLws1bVJ4GOcU/H+mxvAH6gjmeAahGZn8zyWDBxNAJHPI+Pus/lNBFZClwIPAs0qOoJcAIOMDd9JUuJfwf+ARh1H9cCXao64j7Oxc98OdAGfNdt3vuWiJSRw5+1qh4D7gZewwki3cBWcv+zDon12ab8GmfBxCFRnsvpMdMiUg78BPiYqvakuzypJCLXAa2qutX7dJRNc+0zLwAuAr6pqhcCfeRQk1Y0bh/BDcAyYAFQhtPEM16ufdaTSfn5bsHEcRRY5Hm8EDieprKknIgU4gSSH6nqT92nT4Wqve7/rekqXwpcAVwvIodwmjCvwqmpVLtNIZCbn/lR4KiqPus+fgAnuOTyZ30N8KqqtqnqMPBT4A3k/mcdEuuzTfk1zoKJ43lglTviw4fTYfdQmsuUEm5fwbeBl1X1i55fPQS8z/35fcCDM122VFHV21R1oaouxflsn1DVPwK2AO90N8upfQZQ1ZPAERFZ4z51NbCbHP6scZq3LheRUvdcD+1zTn/WHrE+24eAP3VHdV0OdIeaw5LFZsC7ROTtOHer+cB3VPWzaS5SSojIG4HfADuI9B/cjtNvcj+wGOcL+YeqOr5zL+uJSBPwf1X1OhFZjlNTmQO8CPyxqg6ls3zJJiIX4Aw68AEHgQ/g3ETm7GctIp8C3o0zcvFF4EM4/QM59VmLyL1AE06q+VPAJ4GfE+WzdQPr13BGf/UDH1DVlqSWx4KJMcaYRFkzlzHGmIRZMDHGGJMwCybGGGMSVjD5JtmlurpaV65cme5iZIS+vj7KysrSXYyMYMciwo5FhB2LiK1bt55OZA34tAYTEfkOEJpQdm6U3wtObqG344xAeH8oFUgsDQ0NtLQkdZBC1mpubqapqSndxcgIdiwi7FhE2LGIEJHDifx9upu5vseZicq8Up6czBhjTOLSGkxiJCrzSnlyMmOMMYnL9D6TWMnJxszcFJGbcWougFN1NeD3++1YuOxYRNixiLBjkTyZHkziSk6mqvcA9wCIiFobqMPagyPsWETYsYiwY5E86e4zmcysSsBojDHZKtODScqTkxljjElcuocGhxOVichRnERlhQCq+h/AIzjDgvfjJidLT0mNMcZMJK3BRFVvmuT3Cnx4hopjjDFmmjK9mcsYY0wWsGBijDEmYRZMjDHGJCyuPhMRaQSWeLd3Z68bY4wxkwcTEfk8zhKYu4Gg+7QCFkyMMcYA8dVMfh9Yk+3rJRtjjEmdePpMDuLO/TDGGGOiiadm0g9sE5HHgXDtRFVvSVmpjDHGZJV4gslD7j9jjDEmqkmDiap+X0R8wGr3qb2qOpzaYhljjMkm8YzmagK+DxzCSQm/SETeZ0ODjTHGhMTTzPVvwFtUdS+AiKwG7gUuTmXBjDHGZI94RnMVhgIJgKq+go3uMsYY4xFPzaRFRL4N/NB9/EfA1tQVyRhjTLaJJ5j8FU4a+Ftw+kyeBL6RykIZY4zJLvGM5hoCvuj+ywkiztLyznIpxhhjEhUzmIjI/ar6LhHZgZOLawxVPT+lJUsCCxrGGDMzJqqZfNT9/7qZKEiqWWAxxpjUiTmaS1VPuD/+taoe9v4D/npmijczRCQcbIwxxkxdPEODN0R57tpkF8QYY0z2mqjP5K9waiDLRWS751cVwNOpLli6WHOYMcZM3UQ1k/8Gfg8nyePvef5drKp/nIw3F5G3icheEdkvIrdG+f37RaRNRLa5/z6UjPc1xhiTXDFrJqraDXQDNwGIyFygGCgXkXJVfS2RNxaRfODrOM1oR4HnReQhVd09btMfq+pHEnmvZPHWWqwGY4wxEZP2mYjI74nIPuBV4Nc4CR9/mYT3vhTYr6oHVTUA3AfckITXTSrrnDfGmMnFMwP+M8DlwGOqeqGIrMetrSSoETjieXwUuCzKdv9HRK4EXgE+rqpHxm8gIjcDN4ceNzc3j/m993Gsn6f7u/HbZRK/35/R5ZtJdiwi7FhE2LFIHpmsmUZEWlR1nYi8BFyoqqMi8pyqXprQG4v8IfBWVf2Q+/hPgEtV9W8829QCflUdEpG/BN6lqldN8roa2qdYzVLjm6im87tsaOZqbm6mqakp3cXICHYsIuxYRNixiBCRraq6brp/H0/NpEtEynFycv1IRFqBkem+ocdRYJHn8ULguHcDVW33PPwv4PNJeF9jjDFJFs88kxuAAeDjwCbgAM6orkQ9D6wSkWXuSo7vYdzywCIy3/PweuDlJLyvMcaYJIsn0WOf5+H3k/XGqjoiIh8BHgXyge+o6i4R+TTQoqoPAbeIyPU4NaEO4P3Jen9jjDHJM9GkxV7GJngU97EAqqqVib65qj4CPDLuuTs9P98G3Jbo+xhjjEmtieaZVMxkQYwxxmSveDrgEZE3AqtU9bsiUgdUqOqrqS3a9KwGcEdnbAk92dQU8+eJtpvKa2SiC7q6oLo63cXICHYsIuxYRNixSJ54Ji1+EvhHIs1NPuD/pbJQxhhjsks880y2ARcCL6jqhe5z2zN1cSybZxJhY+gj7FhE2LGIsGMRkeg8k3iGBgfcq7O6b1g23TebjSwdizFmNognmNwvIv8JVIvInwOP4UwgNMYYY4D45pncLSIbgB5gDXCnqm5OecmyTLY1gRljTDJNGEzcNPGPquo1gAUQY4wxUU3YzKWqQaBfRKpmqDzGGGOyUDzzTAaBHSKyGQinVlHVW1JWKmOMMVklnmDyC/efMcYYE1U8fSYbkrXmuzHGmNwUT59JvZsi3iSBzTsxxuSieJq5DgFPi8hDjO0z+WKqCmWMMSa7xBNMjrv/8gDLJGyMMeYM8Uxa/BSAiFQ4D9Wf8lIZY4zJKvFkDT5XRF4EdgK7RGSriJyT+qLlPus/Mcbkinhyc90D/K2qLlHVJcDfYbm5jDHGeMQTTMpUNbwWlKo2A5Y5OMnG11Ks1mKMySbxdMAfFJE7gB+6j/8YyMhVFo3JNJYA1MwW8dRM/gyoB37q/qsDPpDKQhljjMku8Yzm6gRSkodLRN4GfBnIB76lqneN+30R8APgYqAdeLeqHkpFWYwxxkxfPKO5NotItedxjYg8mugbu6lavg5cC6wFbhKRteM2+yDQqaorgS8Bn0/0fY0xxiRfPH0mdaraFXqgqp0iMjcJ730psF9VDwKIyH3ADcBuzzY3AP/k/vwA8DXxLvIew/iO6/Ed21PdLhmvkYztTPaa6DM2JhfE02cyKiKLQw9EZAnuevAJagSOeB4fdZ+Luo2qjgDdQO34FxKRm0WkRURaklAuY4wxUxRPzeQTwFMi8mv38ZXAzUl472i3Z+ODVDzboKr34MyHYc2aNbp3797ES5dBJhoRFOt36douE8uUin30yrZRWs3NzTQ1NaW7GBnBjkVEojXmeDrgN4nIRcDl7lMfV9XTCb2r4yiwyPN4IU4OsGjbHBWRAqAK6EjCe2eVbLlIGWNmr3iauQDeADS5/y6fcMv4PQ+sEpFlbor79wAPjdvmIeB97s/vBJ6YrL/EmJmmqhbwzaw3ac1ERO4CLgF+5D71URG5QlVvS+SNVXVERD4CPIozNPg7qrpLRD4NtKjqQ8C3gR+KyH6cGsl7EnlPY4wxqRFPn8nbgQtUdRRARL4PvAgkFEwAVPUR4JFxz93p+XkQ+MNE38cYY0xqxdvMVe35uSoVBTHGGJO94qmZ/AvwoohswRlddSVwe0pLZWIa3zZvbfXGmEwQz2iue0WkGaffRIB/VNWTqS6YmToLLImZKFDbsTVmYvGkU3lcVU+o6kOq+qCqnhSRx2eicMYYY7JDzJqJiBQDpUCdiNQQmUBYCSyYgbIZY4zJEhM1c/0F8DGcwPGC5/kenASNxhhjDDBBMFHVLwNfFpG/UdWvzmCZTBKoKs3NzekuRsazvhBjkiOe0VzdIvKn459U1R+koDzGGGOyUDzB5BLPz8XA1TjNXhZMTFay2ogxyRfP0OC/8T4WkSoi68EbY4wxcdVMxusHVie7ICa17G7cGJNK8SR6/F8ia4jkA2cD96eyUMZMR6yAaYMRjEm9eGomd3t+HsGZb3JTaopjZlqu1lhydb+MyVTx9Jn8WkQuAN4LvAt4FfhJqgtmUscutMaYZJtoBvxqnPVDbgLagR8DoqrrZ6hsxhhjssRENZM9wG+A31PV/QAi8vEZKZUxxpisMlGix/8DnAS2iMh/icjVRPJzGWOMMWExg4mq/kxV3w2cBTQDHwcaROSbIvKWGSqfMcaYLDBpCnpV7VPVH6nqdcBCYBtwa8pLZoxLVaMOGoj1vDFm5sW7bC8Aqtqhqv+pqlelqkAmvTL9wp0p5TDGjDWlYJIsIjJHRDaLyD73/5oY2wVFZJv776GZLqeJLZUXdQsYxmSftAQTnGayx1V1FfA4sZvNBlT1Avff9TNXPAPxX9Tt4m+MSVcwuQH4vvvz94HfT1M5jDHGJEG6gkmDqp4AcP+fG2O7YhFpEZFnRMQCTpbw1lSsdmPM7CCp+gKLyGPAvCi/+gTwfVWt9mzbqapn9JuIyAJVPS4iy4EngKtV9UCU7W4Gbgaor6+/+P77LQ8lgN/vp7y8PN3FyAh2LCLsWETYsYhYv379VlVdN92/T1kwmfBNRfYCTap6QkTmA82qumaSv/ke8LCqPjDRdmvWrNG9e/cmr7BZrLm5maampnQXIyPYsYiwYxFhxyJCRBIKJulq5noIeJ/78/uAB8dvICI1IlLk/lwHXAHsnrESGmOMiVu6gsldwAYR2QdscB8jIutE5FvuNmcDLSLyErAFuEtVLZgYY0wGms5KiwlT1XacteTHP98CfMj9+bfAeTNcNGOMMdOQlj6TVBKRXsA6TRx1wOl0FyJD2LGIsGMRYcciYo2qVkz3j9NSM0mxvYl0IuUSEWmxY+GwYxFhxyLCjkWEiLQk8vfp6jMxxhiTQyyYGGOMSVguBpN70l2ADGLHIsKORYQdiwg7FhEJHYuc64A3xhgz83KxZmKMMWaG5VQwEZG3icheEdkvIrNqNUgRWSQiW0TkZRHZJSIfdZ+Pa+2YXCQi+SLyoog87D5eJiLPusfixyLiS3cZZ4KIVIvIAyKyxz0/Xj9bzwsR+bj7/dgpIveKSPFsOS9E5Dsi0ioiOz3PRT0PxPEV91q6XUQumuz1cyaYiEg+8HXgWmAtcJOIrE1vqWbUCPB3qno2cDnwYXf/4107Jhd9FHjZ8/jzwJfcY9EJfDAtpZp5XwY2qepZwOtwjsmsOy9EpBG4BVinqucC+cB7mD3nxfeAt417LtZ5cC2wyv13M/DNyV48Z4IJcCmwX1UPqmoAuA9n3ZRZQVVPqOoL7s+9OBeMRmbp2jEishB4B/At97EAVwGhRKGz4liISCVwJfBtAFUNqGoXs/S8wJlbVyIiBUApcIJZcl6o6pNAx7inY50HNwA/UMczQLWblDemXAomjcARz+Oj7nOzjogsBS4EniX+tWNyzb8D/wCMuo9rgS5VHXEfz5bzYznQBnzXbfL7loiUMQvPC1U9BtwNvIYTRLqBrczO8yIk1nkw5etpLgUTifLcrBuqJiLlwE+Aj6lqT7rLkw4ich3QqqpbvU9H2XQ2nB8FwEXAN1X1QqCPWdCkFY3bH3ADsAxYAJThNOeMNxvOi8lM+fuSS8HkKLDI83ghcDxNZUkLESnECSQ/UtWfuk+fClVP3f9b01W+GXQFcL2IHMJp7rwKp6ZS7TZvwOw5P44CR1X1WffxAzjBZTaeF9cAr6pqm6oOAz8F3sDsPC9CYp0HU76e5lIweR5Y5Y7M8OF0rD2U5jLNGLdP4NvAy6r6Rc+vJl07Jteo6m2qulBVl+KcB0+o6h/hLGXwTnez2XIsTgJHRCS0+NzVOOsCzbrzAqd563IRKXW/L6FjMevOC49Y58FDwJ+6o7ouB7pDzWGx5NSkRRF5O84daD7wHVX9bJqLNGNE5I3Ab4AdRPoJbsfpN7kfWIzzZfpDVR3fCZezRKQJ+L+qep27/PN9wBzgReCPVXUoneWbCSJyAc5ABB9wEPgAzo3krDsvRORTwLtxRj++iLPkRSOz4LwQkXuBJpxMyaeATwI/J8p54Abbr+GM/uoHPuAuERL79XMpmBhjjEmPXGrmMsYYkyYWTIwxxiTMgokxxpiEWTAxxhiTMAsmxhhjEmbBxJgkEhF/ustgTDpYMDHGGJMwCybGpJiI/J67XsaLIvKYiDS4z9e7a0i8ICL/KSKHRaQu3eU1ZjosmBiTek8Bl7uJFu/DyWYMzgzkJ1T1IuBnOLOQjclKBZNvYoxJ0ELgx24iPR/wqvv8G4EbAVR1k4h0pql8xiTMaibGpN5Xga+p6nnAXwDF7vPR0nwbk5UsmBiTelXAMffn93mefwp4F4CIvAWYFeuwm9xkiR6NSSIRGWXsug9fBA4AX8IJKM8Al6hqk4jMBe7FCSK/xslmuywXM9aa3GfBxJg0EZEiIKiqI7JI4qcAAABWSURBVCLyepzVEC9Id7mMmQ7rgDcmfRYD94tIHhAA/jzN5TFm2qxmYowxJmHWAW+MMSZhFkyMMcYkzIKJMcaYhFkwMcYYkzALJsYYYxJmwcQYY0zC/j+RnrzPzm3QpgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1e2741809b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "x = 2*np.sin(np.arange(100)/10.0) # periodic signal\n",
    "x += np.random.randn(len(x)) # corrupted with noise\n",
    "\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot(x, '-s')\n",
    "plt.ylabel('x', fontsize=20)\n",
    "plt.grid(True)\n",
    "plt.xlabel('Time')\n",
    "\n",
    "plt.subplot(2,1,2)\n",
    "c = plt.acorr(x, usevlines=True, normed=True, maxlags=50, lw=2)\n",
    "plt.grid(True)\n",
    "plt.axhline(0, color='black', lw=2)\n",
    "plt.axhline(1/np.exp(1), color='red')\n",
    "plt.ylabel('Autocorrelation')\n",
    "plt.xlim(xmin=0,xmax=100)\n",
    "plt.xlabel('Lag')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>图5.18:显示使用正弦函数生成的100个随机数的图，以及该系列的自相关。</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为了精确地确定相关长度，我们将拟合lag和相关长度之间的插值函数，然后确定自相关变为$e^{-1}$的lag。`plt.acorr`在这些lag处返回滞后和自相关。首先我们将各自为lag和自相关性赋单独的变量。我们也可以打印lag，看看里面是什么。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-50 -49 -48 -47 -46 -45 -44 -43 -42 -41 -40 -39 -38 -37 -36 -35 -34 -33\n",
      " -32 -31 -30 -29 -28 -27 -26 -25 -24 -23 -22 -21 -20 -19 -18 -17 -16 -15\n",
      " -14 -13 -12 -11 -10  -9  -8  -7  -6  -5  -4  -3  -2  -1   0   1   2   3\n",
      "   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21\n",
      "  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39\n",
      "  40  41  42  43  44  45  46  47  48  49  50]\n"
     ]
    }
   ],
   "source": [
    "lags = c[0] # lags\n",
    "auto_corr = c[1] # autocorrelation\n",
    "print(lags)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`acorr`提供正、负lags。我们并不都需要，我们可以通过给出一个布尔格式的索引来解决这个问题，并通过使用`lags>=0`语句获得索引。我们还删除了对应于负lags的自相关数组。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "auto_corr = auto_corr[lags>=0]\n",
    "lags = lags[lags>=0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在，我们需要两个点的自相关，一个只在阈值以上，另一个个在阈值以下。我们通过计算自相关超过阈值的次数来得到它们的索引。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9\n"
     ]
    }
   ],
   "source": [
    "n = sum(auto_corr>np.exp(-1))\n",
    "print(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "一个点是在第八个索引处，和另一个点在第九个索引处。现在，当x的自相关等于阈值时，我们可以用`interp1d`来求精确的lag值。这里提供的相关长度为8.74。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(8.742091971219681)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.interpolate import interp1d\n",
    "f = interp1d([auto_corr[n], auto_corr[n-1]], [lags[n], lags[n-1]])\n",
    "corr_len = f(np.exp(-1))\n",
    "corr_len"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
